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ABSTRACT 


A  nonlinear  elastic  (hypoelastic)  constitutive  relation  is 
proposed  for  analysis  of  plane  and  axisymmetric  reinforced  and/or 
prestressed  concrete  structures.  The  constitutive  relation  is  based 
on  the  equivalent  uniaxial  strain  concept.  A  now  characterization 
of  Poisson's  ratio  is  introduced.  Strain  softening  behavior  is 
assumed  in  tension  and  post  failure  conditions  are  imposed  on  three 
dimensional  ultimate  strength  and  corresponding  equivalent  uniaxial 
strain  surfaces.  Special  isoparametric  finite  elements  are  developed 
for  representation  of  reinforcing  bars  and  prestressing  layers. 

The  finite  element  model  as  well  as  the  proposed  constitutive 
relations  are  incorporated  in  a  finite  element  program  (FEPARCS5) 
for  nonlinear  analysis  of  axisymmetric  or  plane  reinforced  and/or 
prestressed  concrete  structures.  Program  FEPARCS5  is  then  used  to 
analyze  a  finite  element  model  of  an  axisymmetric  prestressed 
concrete  test  structure  resembling  a  secondary  containment  building 
under  increasing  internal  pressure.  The  results  of  the  analysis 
are  then  compared  to  the  results  of  the  test  structure  and  a 
theoretical  elastic  plastic  analysis  of  the  same  test  structure. 
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CHAPTER  ONE 


INTRODUCTION 

1 . 1  Background  to  the  Problem 

In  the  field  of  nuclear  energy,  extreme  risks  are  associated 
with  accident  conditions  such  as  overpressure,  earthquake  loads  or 
catastrophic  impact.  Therefore,  the  ability  to  predict  the  entire 
response  of  a  structure  such  as  a  secondary  containment  building 
or  a  primary  containment  vessel  and  to  identify  the  conditions  under 
which  limit  states  occur  in  order  to  assess  possible  damage  is 
highly  desirable.  This  task  calls  for  complicated  nonlinear 
analyses  which  are  generally  expensive.  Often  economic  constraints 
lead  the  analyst  either  to  abandon  or  to  oversimplify  the  problem. 
These  considerations  call  for  the  development  of  more  efficient  and 
versatile  solution  procedures  (Almroth,  Stren  and  Brogan,  1979) . 

There  are  several  aspects  to  the  nonlinear  analysis  of 
such  structural  problems.  Some  of  the  major  aspects  are  the  con¬ 
stitutive  modelling  of  the  material,  the  finite  element  modelling 
of  the  structure  and  the  loads,  and  the  numerical  solution  tech¬ 
nique  of  the  nonlinear  problem.  In  addition,  it  is  important  to 
realize  that  there  is  a  strong  interaction  between  these  aspects 
(Bathe  and  Ramaswamy,  1979) . 

Since  1974  a  research  program  has  been  underway  at  the 
University  of  Alberta,  Edmonton  sponsored  by  the  Atomic  Energy 
Control  Board  of  Canada  to  investigate  the  effect  of  over-pressure 
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on  Gentilly-2  type  secondary  containment  structures  which  house 
CANDU  nuclear  reactors.  An  advanced  elastic  plastic  constitutive 
relation  for  biaxial  behavior  of  concrete  has  been  developed 
(Epstein  and  Murray,  1978  and  Murray,  et  al.,  1978)  and  implemented 
by  Murray,  Chitunyanondh  and  Wong  (1978)  in  a  modification  of  the 
B0S0R5  code  (Bushnell,  1973).  Subsequently,  a  series  of  tests 
were  conducted  on  reinforced  and  prestressed  concrete  wall  segments 
under  biaxial  and  uniaxial  tension  to  assess  the  performance  and 
parameters  of  the  constitutive  relation.  In  addition,  a  reinforced 
and  prestressed  test  structure  composed  of  a  cylinder  and  a  dome 
was  built  and  tested  under  internal  pressure.  The  test  results 
have  been  compared  to  the  modified  B0S0R5  analysis  of  the  same 
test  structure  in  order  to  assess  the  capability  of  the  program 
and  the  constitutive  relation  to  predict  the  behavior  of  the  test 
structure . 

Although  this  thesis  is  not,  formally,  a  part  of  that 
research  program,  the  motivation  behind  the  work  reported  herein,  is 
to  develop  a  parallel  sophisticated  capability  founded  upon  alterna¬ 
tive  technology.  A  three  dimensional  constitutive  relation  is  deve¬ 
loped  to  take  account  of  thick  and  thin  shell  action  and  an  axisym- 
metric  finite  element  thickshell  model  is  used  to  accommodate  it.  Both 
are  implemented  in  a  nonlinear  finite  element  program  for  analysis 
of  axisymmetric  reinforced  and/or  prestressed  concrete  structures. 
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1.2  Scope  and  Objectives  of  Thesis 

The  scope  of  this  thesis  is  the  nonlinear  static  analysis 
of  three  dimensional  (axisymmetric)  reinforced  and/or  prestressed 
structures.  Displacements  are  small,  rotations  are  negligible  and 
strains  are  assumed  to  be  infinitesimal.  High  temperature  and  creep 
effects  are  outside  the  scope  of  this  work. 

The  objectives  of  the  study  can  be  listed  as 

1.  To  develop  a  three  dimensional  nonlinear  elastic 
constitutive  relation  for  concrete. 

2.  To  formulate  a  finite  element  model  capable  of 
representing  axisymmetric  behavior  of  reinforced 
and/or  prestressed  concrete  structures. 

3.  To  develop  a  nonlinear  finite  element  program  for 
analysis  of  such  structures. 

4.  To  demonstrate  the  capabilities  of  this  program 
through  application  to  a  complicated  structure  such 
as  a  prestressed  concrete  containment  structure  and 
to  identify  the  limit  states  associated  with  over¬ 
pressure  loading  of  this  structure. 

1 . 3  Organization  of  Thesis 

Chapter  Two  contains  a  literature  review  of  some  of  the 
existing  constitutive  relations  for  multiaxial  behavior  of  concrete, 
and  the  development  and  description  of  a  proposed  nonlinear  elastic 
constitutive  relation  for  three  dimensional  (axisymmetric)  behavior 


of  concrete. 
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In  Chapter  Three  a  finite  element  model  for  incremental 
displacement  analysis  of  axisymmetric  reinforced  and/or  prestressed 
concrete  structures  is  presented.  The  necessary  incremental  vari¬ 
ational  principles  are  discussed.  Finite  elements  representing 
concrete,  meridional  reinforcing  and  circumferential  reinforcing 
are  formulated  and  the  associated  boundary  conditions  and  the  work 
equivalent  loads  are  derived. 

A  finite  element  program  (FEPARCS5)  for  nonlinear  analysis 
of  axisymmetric  reinforced  and/or  prestressed  concrete  structures 
is  described  in  Chapter  Four.  The  solution  techniques  incorporated 
in  the  program,  as  well  as  special  capabilities  such  as  initial 
stress  and  post-tensioning  simulation,  are  discussed  and  the  flow  of 
operations  is  presented. 

A  preliminary  investigation  of  the  capabilities  of  program 
FEPARCS5  is  carried  out  in  Chapter  Five.  Chapter  Six  contains  an 
analysis,  using  program  FEPARCS5,  of  a  prestressed  containment 
structure  under  internal  pressure  which  was  built  and  tested  to 
failure,  in  the  I.F.  Morrison  Structural  Laboratory.  The  test 
structure  and  procedure  are  briefly  described.  The  finite  element 
model  of  the  test  structure  is  presented.  The  results  of  the 
analysis  are  compared  with  the  test  results  as  well  as  with  the 
results  of  an  elastic  plastic  analysis.  Finally  some  of  the  limit 
states  associated  with  the  loading  program  are  identified. 

In  Chapter  Seven  conclusions  are  drawn  on  the  performance 
of  the  constitutive  model  and  the  capabilities  of  program  FEPARCS5. 


CHAPTER  TWO 


CONSTITUTIVE  THEORY 


2 . 1  Introduction 

A  multiaxial  constitutive  relationship  is  an  essential  com¬ 
ponent  in  any  finite  element  nonlinear  analysis  of  reinforced 
concrete  structures.  Unfortunately,  the  behavior  of  concrete  under 
multiaxial  states  of  stress  is  complex  both  in  the  strength  and  in 
the  deformation  domains.  While  information  on  uniaxial  response  of 
concrete  is  abundant,  biaxial  and  triaxial  responses  are  not  yet 
fully  understood.  This  situation  is  not  improved  by  the  scarcity 
of  reliable  data  on  which  to  base  analytical  models.  There  is  a 
general  lack  of  strain  data  particularly  near  and  beyond  peak 
strengths.  Cracking  and  post-crushing  behavior  are  major  problems 
since  a  large  part  of  the  response  to  failure  of  a  reinforced  con¬ 
crete  structure  must  of  necessity  be  traced  after  part  or  most  of 
the  structure  have  cracked  and  some  concrete  may  have  crushed. 

The  scope  of  analysis  of  concrete  structures  can  range 
from  linear  small  displacement  infinitesimal  strain  analysis  to 
large  displacement  finite  strain  analysis.  Therefore,  choice  of 
constitutive  relations  suitable  to  the  type  of  analysis  enhances 
efficiency.  However,  the  physical  nonlinearity  of  concrete  is 
always  dominant.  Approaches  to  the  material  model  range  from 
elastic-plastic  to  nonlinear  elastic.  One  popular  approach  to 
modelling  the  behavior  of  concrete  under  multiaxial  states  of 
stress  is  to  consider  the  material  to  be  orthotropic  nonlinear 
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elastic  and  to  organize  the  constitutive  model  around  some  equiva¬ 
lent  uniaxial  relation.  Since  1972  a  number  of  studies  have  been 
published  along  that  line;  Liu,  Nilson  and  Slate  (1972),  Coon  and 
Evans  (1972),  Kupfer  and  Gerstle  (1973),  Romstadt,  et  al.  (1974), 
Darwin  and  Pecknold  (1974,  1977a  and  1977b),  Same  (1974),  Link 
(1976)  and  Bashur  and  Darwin  (1978) . 

In  this  chapter  a  brief  review  of  some  of  these  models  is 
presented  and  a  nonlinear  three  dimensional  (axisymmetric)  consti¬ 
tutive  relation  for  concrete  is  developed. 

2 . 2  Literature  Review 

Truesdell  (1955)  defines  hypoelastic  materials  as  those  for 
which  the  rate  of  stress  is  a  function  of  the  rate  of  deformation 
and  the  stress  history.  Argyris,  et  al.  (1976)  and  Schnobrich  (1977) 
classify  the  models  referred  to  in  Section  2.1  as  hypoelastic.  The 
following  discussion  is  confined  to  this  approach  since  the  model 
proposed  by  the  writer  falls  under  the  same  classification.  The 
underlying  assumption  of  these  models  besides  the  definition  of 
Truesdell  is  that  the  material  is  nonlinear  elastic.  Therefore, 
tangent  or  secant  constitutive  equations  are  used  in  the  form  of 
the  generalized  Hooke's  Law.  The  problem  then  is  to  determine  the 
variation  of  the  moduli  involved  in  a  given  form  of  Hooke's  Law 
throughout  the  loading  history.  In  the  absence  of  time  dependent 
effects  and  high  temperatures,  these  moduli  become  functions  of 
the  level,  type  and  ratios  of  the  stresses  and/or  the  strains. 
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Liu,  Nilson  and  Slate  (1972)  have  proposed  a  biaxial  ortho¬ 
tropic  model  of  concrete.  In  this  model  a  stress-strain  curve  which 
takes  into  consideration  Poisson's  ratio  and  the  ratio  of  principal 
stresses  is  defined  as 

a.  =  e.E  /{  (1  -  va)  [1+  (e/e  .  (1  -  va)  -  2)  (e ./e  . ) 

1  1  O  SI  1  Cl 

+  (£./e  .)2]}  ,  i  =  1/2  (no  sum)  (2.1) 

1  ci 


where,  E  is  the  secant  modulus  at  peak  strength,  £  is  the  strain 
s  c 

at  peak  strength  and  a  is  the  principal  stress  ratio.  The  incre¬ 
mental  constitutive  equation  is  written  as 


pH 

<3 

A  Elb/E2b  Av  0 
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Aq2 

= 

Av  A  o 

Ae2 

At  12 

°  °  (Elb  E2b)/(E‘b  +  E2b+2E2bJ) 

CM 

r-H 

>- 

<3 

(2.2) 


where , 


A 


Elb/(Elb/E2h  " 


2b 


(2.3) 


and  Ej  and  E2^  are  effective  tangent  moduli,  derived  by  differen¬ 
tiation  of  Eq.  2.1  with  respect  to  the  strain  £_^.  Eq.  2.1  is  used 
in  compression  only,  while  in  tension  a  straight  line  relation 
ending  with  a  tension  cut  off  criterion  is  used. 

Kupfer  and  Gerstle  (1973)  have  suggested  a  secant  isotropic 
model  for  biaxial  behavior  of  concrete  in  which  the  shear  and  bulk 

moduli,  respectively  G  and  K  ,  are  used  as  expressed  in  the 

s  s 


equation. 
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(3K  +  G  ) 
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0 

£2 
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0 

0 

(3K  +  4G  )/4 
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(2.4) 


Kupfer  and  Gerstle  suggest  relationships  for  the  shear  and  bulk  moduli 
as  functions  of  the  octahedral  shear  strain  as 


G  =  G  (1  -  a  (T  /f  )m) 
s  o  oc  cu 


(2.5a) 


K  =  K  (G  /G  )  exp  (cyP  ) 
s  o  s  o  oc 


(2.5b) 


where  G^  and  Kq  are  the  initial  shear  and  bulk  moduli  respectively, 

T  and  Y  are  the  octahedral  shear  stress  and  strain,  f  is  the 
oc  oc  cu 

uniaxial  compressive  strength  and  a,  m,  c  and  p  are  constants.  These 

authors  have  used  the  failure  surface  proposed  by  Kupfer,  Hilsdorf 

and  Riisch  (1969)  as  a  biaxial  failure  criterion.  This  type  of 

analysis  agrees  well  with  test  data  in  the  compression 

zone  but  is  not  successful  in  other  zones  (Darwin  and  Pecknold, 

1977b  and  Kupfer  and  Gerstle,  1973). 

Cedolin ,  Crutzen  and  Poli  (1976  and  1977)  have  suggested 

an  extension  of  the  same  concept  to  three  dimensional  cases.  However, 

the  bulk  modulus  in  this  case  is  assumed  to  be  a  function  of  the 

normal  octahedral  strain  £  .  Their  recommended  expressions  are 

o 


G  =  G  (.81  (2  Toc/*002)  _ 2y  +.9) 
so  oc 


(2.6a) 


K  =  K  (.85  (2.5  e°/-0014)  +  .is) 
s  o 


(2.6b) 
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While  Argyris,  et  al.  (1974)  admit  that  there  may  be  a 
relation  between  the  deviatoric  and  hydrostatic  moduli,  they 
state  that  the  assumption  of  isotropy  is  limiting  under  general 
conditions  such  as  nonproportional  load  paths  and  cyclic  loading. 

Romstadt,  et  al.  (1974)  has  defined  a  strain  space  which  is 
divided  into  several  progressive  damage  zones.  As  the  material 
strains  from  one  damage  zone  to  another.  Young's  modulus  and 
Poisson's  ratios  assume  different  values.  The  variation  of  Young's 
modulus  takes  into  account  the  previous  history  of  damage.  Although 
this  concept  is  very  interesting,  it  has  not  been  pursued  by  other 
investigators . 

Link  (1976  and  1977)  has  developed  a  number  of  functions  to 
describe  Young's  moduli  and  Poisson's  ratios  to  be  used  in  a  two 
dimensional  orthotropic  material  model.  These  functions  employ  the 
level  and  ratio  of  principal  stresses  as  independent  variables. 

The  functions  are  extremely  complicated  and  have  not  been  extended 
to  three  dimensional  behavior. 

Geistfeldt  (1977a  and  1977b)  has  suggested  a  hyperelastic 
approach  to  obtain  the  required  elastic  moduli  for  an  isotropic 
biaxial  behavior  of  concrete.  Geistfeldt  (1977b)  suggests  that 
limited  extension  to  three  dimensions  can  be  achieved. 

Darwin  and  Pecknold  (1974,  1977a  and  1977b)  have  developed 
the  concept  of  equivalent  uniaxial  strains  for  orthotropic  biaxial 
behavior  of  concrete.  The  incremental  constitutive  relation  is  of 


the  form 


♦ 
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Aai 

Ej  v/e  iE2 

0 

Aei 

Aa2 

=  1/(1 -v2) 

v/EiE2  e2 

0 

As2 

Ax  12 

0  0 

Gi2(l  -V2) 

Ay  1 2 

(2.7) 

The  equivalent  uniaxial  strain  increment  is  defined  as 


Ae  .  =  Aa./E.  (no  sum)  ,  i  =  1,2  (2.8) 

Ul  1  1 

Eqs.  2.7  and  2.8  are  defined  in  the  principal  axes  of  orthotropy 
which  are  assumed  to  lie  along  the  axes  of  principal  stresses.  The 
accumulation  of  the  increments  of  equivalent  uniaxial  strains  yields 
the  total  equivalent  uniaxial  strains  which  generally  do  not  form  a 
second  order  tensor  and,  therefore,  are  not  transformable. 


£  .  =  ZAa . /E .  (no  sum)  ,  i  =  1,2  (over  loadpath)  (2.9) 

Ul  11 

Darwin  and  Pecknold  have  chosen  the  uniaxial  compressive 
stress  strain  relation  of  Saenz  (1964)  to  act  as  an  instantaneous 
stress  equivalent  uniaxial  strain  relation  in  compression.  A 
straight  line  relation  ending  in  a  cut  off  has  been  used  in  tension. 

The  parameters  of  peak  strength  appearing  in  the  relation  are 

•  • 

obtained  from  a  modified  form  of  the  Kupfer,  Hilsdorf  and  Rusch 
failure  envelope  (1969) . 

It  is  noted  that  most  investigators  have  used  a  tension 
cutoff  criterion.  The  proponents  of  the  gradual  softening  approach 
in  tension  claim  that  cracking  (when  treated  in  a  smearing  fashion) 
cannot  occur  over  every  point  in  a  region  where  stiffness  is  rep¬ 
resented  by  an  integral.  Therefore,  the  material  must  retain  a 
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a  certain  measure  of  stiffness  represented  by-  a  descending  branch 
(Scanlon  and  Murray,  1974,  and  1972,  Murray,  et  al . ,  1978  and 
Elwi  and  Murray,  1979) .  The  descending  branch  in  tension  is  observed 
in  indeterminate  tests,  e.g.  Evans  and  Marathe  (1968),  and  its 
characteristics  must  be  functions  of  the  steel  ratio  in  real 
structures . 

The  above  has  been  a  brief  discussion  of  some  of  the 
proposed  models  for  the  behavior  of  concrete  under  multiaxial 
states  of  stress,  and  has  been  confined  to  hypoelastic  proposals. 

Many  other  contributions  of  equal  importance  have  been  published. 

Coon  and  Evans  (1972)  have  introduced  a  true  hypoelastic  approach; 
Murray  (1979)  has  developed  octahedral  based  tangential  stiffness 
matrices;  Same  (1974)  has  extended  the  hypoelastic  approach  to 
•three  dimensions;  Ottosen  and  Andersen  (1975  and  1977)  have  proposed 
yet  different  models;  Ottosen  (1979)  has  recently  proposed  a  secant 
based  constitutive  relation  which  tackled  the  post  failure  behavior. 
On  the  other  hand,  elastic  plastic  models  have  been  used  and  are 
still  being  developed.  Among  these  are  the  Chen  and  Chen  (1975) 
three  dimensional  elastic  plastic  model,  the  investigation  by 
Mroz  (1972)  of  nonassociated  flow  rules,  the  three  parameter 
elastic  plastic  biaxial  constitutive  theory  suggested  by  Epstein 
and  Murray  (1978)  and  the  further  developments  by  Murray,  et  al. 

(1978) .  In  a  separate  category  the  endochronic  theory  developed  by 


Bazant  and  Coworkers  (1975  and  1976)  must  be  mentioned. 


•  _ 
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2 . 3  The  Proposed  Material  Model 

2.3.1  Introduction 

A  constitutive  relationship  based  on  a  hypoelastic 
orthotropic  approach  is  proposed  to  model  the  behavior  of  a  three 
dimensional  (axisymmetric)  concrete  continuum.  The  model  assumes 
small  displacements,  infinitesimal  strains  and  negligible  rotations. 
No  rate,  temperature  or  creep  effects  are  provided  for.  It  is 
defined  in  the  form  of  an  incremental  stress-strain  constitutive 
equation  in  which  the  material  parameters  are  obtained  from  stress- 
equivalent  uniaxial  strain  relationships.  The  model  draws  heavily 
on  the  earlier  work  of  Darwin  and  Pecknold  (1974,  1977a  and  1977b) , 
of  Saenz  (1964)  and  of  Wiliam  and  Wamke  (1975).  In  the  following, 
the  constitutive  matrix,  the  stress  equivalent  uniaxial  strain 
relation,  the  constitutive  or  material  parameters  and  the  ultimate 
strength  and  corresponding  equivalent  uniaxial  strain  criteria, 
incorporated  into  the  model,  are  discussed. 

2.3.2  The  Constitutive  Equation 

The  proposed  constitutive  relation  is  intended  for  use  in 
the  analysis  of  axisymmetric  structures.  Therefore,  the  constitu¬ 
tive  matrix  is  4  x  4.  When  orthotropy  is  assumed,  the  number  of 
independent  variables  in  the  constitutive  matrix  is  restricted  to 
ten.  When  the  relation  is  referred  to  the  principal  axes  of  ortho¬ 
tropy  it  can  be  written  in  the  form  of  an  incremental  Hooke's 


Law  as  follows 
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It  will  be  assumed  that  the  constitutive  matrix  will  remain 
symmetric  throughout  the  analysis.  This  symmetry  gives  rise  to 
the  following  relations 


V  12  Ei 

=  V2  1 

e2 

(2.11a) 

V  1  3  Ei 

=  ^3  1 

e3 

(2.11b) 

V  2  3  E2 

=  V3  2 

e3 

(2.11c) 

substituting  Eqs .  2.11  into  Eq.  2.10  and  rearranging,  the  following 
symmetric  form  of  Eq.  2.10  can  be  obtained, 
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where 
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(2 . 14d) 

Eqs .  2.13  and  2.14  are  the  basis  of  the  constitutive 
relation.  The  symmetry  assumed  above  is  by  no  means  trivial 
and  will  be  enlarged  upon  later. 


2.3.3  The  Equivalent  Uniaxial  Strain  Concept 

Having  defined  the  form  of  the  incremental  constitutive 
equation  it  remains  to  determine  the  variation  of  the  tangential 
hypoelastic  moduli.  For  this  purpose,  the  concept  of  equivalent 
uniaxial  strains  which  has  been  developed  by  Darwin  and  Pecknold 
(1974,  1977a,  1977b)  is  used  with  slight  modifications.  This 
concept  is  briefly  described  as  follows. 

Let  Eq.  2.13  be  written  as 
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E 1 B 1  3 

0 

de  1 

dG  2 

E2B2  1 

E2B22 

E2B2  3 

0 

de2 

CO 

to 

E3B3  1 

E  3B  3  2 

E  3B  3  3 

0 

de3 

dT  1 2 

0 

0 

0 

G12 

dy  1 2 

in  which  the  coefficients  B__  are  defined  by  identifying  the  matrix 
terms  of  Eq.  2.15  with  the  corresponding  terms  of  Eq.  2.13. 


Carrying  out  the  multiplications  of  Eq.  2.15,  yields 


- 


* 
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dO i  =  Ei  (Biid£i  +  B  i  2  d£  2  +  B i 3  d£  3 ) 

dO 2  =  E2  (B2ld£i  +  B22d£2  +  B23d£3) 

dO 3  =  E3  (B3id£i  +  B  3  2  d£  2  +  B  3  3  d£  3 ) 

dTi2  =  G 1 2  dy  1 2 

Eqs.  2.16  is  written  in  matrix  form  as  follows 


do  1 

dO  2 

dO  3 

dT  1  2 

Ei 

0 

0 

0 


0 

e2 

0 

0 


0 

0 

E3 

0 


0 

0 

0 

Gi  2 


d£ 


d£ 


lu 


2  u 


d£ 

3  U 

dyi2 


(2.16a) 
(2.16b) 
(2.16c) 
( 2 . 16d) 


(2.17) 


where  <  d£^  >  is  defined  as  the  equivalent  uniaxial  strain  increment 
and  is  written  in  terms  of  the  actual  strain  increment  as 

3 

d£  . =  £  b.  .  d£ .  ,  i  =  1,3  (2.18) 

U1  j=i  13  3 

The  equivalent  uniaxial  strain  increment  can  be  evaluated  from 
Eq.  2.17  in  the  simple  form 


d£  .  =  dO./E.  (no  sum)  ,  i  =  1,3  (2.19) 

ui  11 

and  the  total  equivalent  uniaxial  strain  may  be  determined  by 
integrating  Eq.  2.19  over  the  load  path 


£  .  =  dO./E.  (no  sum)  ,  i  =  1,3  (2.20) 

ui  11 

(over  load  path) 


Eqs.  2.16  are  defined  in  the  principal  axes  of  orthotropy. 


Darwin  and  Pecknold  assumed  that  the  principal  axes  of  orthotropy 


. 
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follow  the  principal  axes  of  stresses.  However,  in  axisymmetric 
shell  structures  the  meridional  and  circumferential  stress  directions 
predominate.  Therefore,  the  proposed  model  assumes  that  the  princi¬ 
pal  axes  of  orthotropy  at  a  point  are  fixed  in  the  local  meridional, 
circumferential  and  normal  directions .  The  equivalent  uniaxial 
strains  will  be  treated  in  the  same  manner  as  real  strains  and  will 
be  transformed  using  the  usual  coordinate  transformation  methods.  A 
physical  interpretation  for  the  equivalent  uniaxial  strains  and  the 
diagonal  matrix  of  Eq.  2.17  is  attempted  in  the  following.  The 
increment  equivalent  uniaxial  strain  of  Eq.  2.19  is  the  increment 
of  the  strain  in  direction  i  that  the  material  would  exhibit  if 
subjected  to  a  stress  increment  dG_^  while  all  other  stress  increment 
contributions  are  equal  to  zero.  In  other  words,  the  matrix  of 
Eq.  2.17  is  the  constitutive  matrix  of  a  fictitious  material  with 
zero  Poisson's  ratios. 

2.3.4  The  Equivalent  Uniaxial  Strain-Stress  Relation 

In  the  hypoelastic  theory  developed  by  Truesdell  (1955) 
the  stress  strain  relation  follows  from  the  incremental  constitutive 
equations.  In  the  material  model  proposed  herein  it  is  assumed 
that  the  total  stresses  are  functions  of  the  current  equivalent 
uniaxial  strains. 

For  the  purposes  of  this  work,  the  uniaxial  compressive 
stress  strain  relationship  of  Saenz  (1964)  is  generalized  in  a 


.r  ; 
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manner  similar  to  that  proposed  by  Darwin  and  Pecknold  (1974,  1977a 
and  1977b)  and  used  to  describe  compressive  as  well  as  tensile 
responses  on  the  ascending  part  of  the  stress-equivalent  uniaxial 
strain  curve.  Writing  Saenz's  relationship  in  terms  of  the  equi¬ 
valent  uniaxial  strain  yields 


Eoi£ui//^+  (£ui/£ci^  +  ^eui^£ci^  ^  (2.21) 


a. 

l 


where 


Rg  =  E0i/Egi  (no  sum) 


(2.22a) 


Ck,.:  / £_,.•  (no  sum) 


(2.22b) 


The  type  of  curve  described  by  Eq.  2.21  is  illustrated  in  Fig.  2.1 

on  which  the  variables  of  Eqs .  2.21  and  2.22b  are  shown.  More 

specifically,  is  the  initial  modulus  of  elasticity  in  direction  i 

and  CT  •  and  £  •  are  the  maximum  stress  associated  with  direction  i 
ci  ci 

and  the  corresponding  equivalent  uniaxial  strain  respectively. 

The  descending  branch  of  Eq.  2.21  is  too  steep  for  lightly 
reinforced  concrete  in  tension.  Therefore,  in  this  study  a  bilinear 
descending  branch  is  adopted,  as  follows  (no  sum) 


a. 

l 


aci  a2  (32  -  eui/£ci)  <  £ui/£ci  1  $2  (2.23b) 


a. 

1 


0 


(2.23c) 


where  ai , 


a2 ,  3i  and  32  are  shown  on  Fig.  2.2. 


■ 
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In  compression,  Eq.  2.21  is  retained  up  to  £  .  =  2e  .  at  which 

ui  ci 

point  it  is  assumed  that  concrete  crushes.  The  complete  curves 
for  compression  and  tension  are  shown  on  Figs.  2.1  and  2.2 
respectively . 

Eq.  2.16d  shows  that  the  shear  strain  can  be  accumulated 
directly.  For  the  purposes  of  this  study  the  Saenz  relation  was 
adopted  for  shear  and  is  written  as 

T12  =  Go12Y12/[1  -  (Rg  -  2)  (Yi2/YCl2)  +  (Y12  YCi2 ) 2]  (2.24) 

where 


rg  - 

gs12/go12 

(2.25a) 

G 

Sl  2 

=  t  ,/y 

C  1  2  C  1  2 

(2.25b) 

and  G  . „  is  the  initial  shear  modulus,  and  T  and  Y  are  the 
oi2  c i 2  c i 2 

maximum  shear  stress  and  the  corresponding  shear  strain  respectively. 
Eq.  2.24  is  adopted  up  to  the  peak  of  the  shear  curve.  For  the 
descending  branch  of  the  shear  response  a  straight  line  is  used 
as  follows 


1 2 


=  T 


Cl  2 


(l-a3  (Yi2/YC12  -  D)  1  <  t12/t 


ci  2 


(2.26) 


where  a  3  is  a  constant.  The  shear  stress  strain  curve  is  shown  in 
Fig.  2.3.  In  this  study  shear  deformations  are  expected  to  be  low 
in  the  particular  coordinate  system  chosen.  Therefore,  the  difference 
between  Eq.  2.24  and  the  equations  of  Kupfer  and  Gerstle  (1973) 
or  Cedolin,  et  al_  (1976  and  1977)  is  expected  to  be  of  no  importance 


and  the  point  is  not  further  pursued. 


. 
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2.3.3  Young's  Moduli  and  the  Shear  Modulus 

Eq.  2.21  serves  to  define  the  incremental  elastic  moduli 
of  Eq.  2.17,  for  by  Eq.  2.19 

E_^  =  dO\/d£u^  (no  sum)  (2.27) 


Differentiating  Eq.  2.21  with  respect  to  £  yields 

1U 


Ei  “  Eoi(l-(£ui/£oi)2/[1-(RE-2)(eui/£ci)  +  (0ui/Eoi)2]2 

(no  sum)  (2.28) 


which  defines  Young's  moduli. 

The  incremental  shear  modulus  G12  can  be  obtained  from 
Eq.  2.24  using  similar  reasoning  as 


2  i  2 


G12  =  G012  (1-  (Y.2/TCl2)")/[l-  (Rg-2)(Y12/Tc12)  +  <Y12/Ycl2n 

(2.29) 


Eqs .  2.28  and  2.29  are  confined  to  the  ascending  branches  of  the 
respective  responses,  except  in  compression  where  Eqs.  2.28  is  still 
used  for  the  descending  branch.  In  tension  and  shear  however,  the 
required  moduli  are  assumed  to  be 


E.  = 

1 

-  <*! 

Esi 

1  < 

£  .  /£  .  < 

Ul.  Cl 

01 

(2.30a) 

E  . 

1 

-  02 

Esi 

3i  < 

£  .  /£  .  < 
Ul  Cl 

32 

(2.30b) 

G 12  = 

-  O3 

GS12 

1  < 

y.2/yc12 

(2.30c) 

The  choice  of  ai  ,  O2 ,  $i,  32  and  013  depends  on  the  amount 
of  steel  in  the  neighbourhood  of  the  point  under  consideration.  In 
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order  that  a  tangential  stiffness  approach  be  successful,  it  is 
recommended  that  the  overall  uniaxial  stiffness  of  a  region 
including  steel  be  required  to  be  positive  definite  ( Chi tnuy anon dh, 
et  al. ,  1979) . 

2.3.6  Poisson's  Ratio 

Prior  to  implementing  the  incremental  stress  strain  relation¬ 
ship  it  is  necessary  to  determine  the  values  of  Poisson's  ratios 
appearing  in  Eq.  2.13.  It  is  assumed  herein  that  a  strain  de¬ 
pendent  Poisson's  ratio  may  be  applied  to  each  equivalent  uni¬ 
axial  strain,  that  is,  three  independent  Poisson's  ratios  are 
postulated  in  the  form 

V.  =  V  f  (e  ./£  .)  (no  sum)  (2.31) 

1  o  Ul  Cl 


in  which  V  is  the  initial  Poisson's  ratio.  Eqs .  2.14a  to  2.14c 
o 

may  now  be  written  as 


=  V!  V2 


1 2  3  =  V2  V3 


rt  =  v  v 

1  3  1  3  1 


(2.32a) 


(2.32b) 

(2.32c) 


In  compression,  the  function  fCe^/e^.  )  appearing  in 
Eq.  2.31  has  been  determined  from  the  uniaxial  compression  data  of 
Kupfer,  Hilsdorf  and  Rusch  (1969)  by  a  least  squares  fit  of  a  cubic 
polynomial.  This  results  in  the  approximation 

V  =  V  (1.0000  +  1.3763  (e/e  )  -  5.3600  (e/e  )2 
o  cu  cu 

+  8.5860  (e/e  )  3) 
cu 


(2.33) 


' 
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in  which  £  is  the  strain  in  the  direction  of  loading  and  £  is  the 

cu 

strain  at  ultimate  strength.  Substituting  £  .  and  £  .  for  £  and 

ui  ci 

£cu  respectively,  Eq.  2.33  assumes  the  form  of  Eq.  2.31. 

Micro-cracks  exist  in  concrete  along  the  aggregate  mortar 
interfaces  even  at  the  virgin  state  (Hsu,  et  al . ,  1963).  At  30%  of 
ultimate  strength  these  cracks  increase  in  number,  width,  and  length. 
At  approximately  75%  of  ultimate  strength  these  cracks  penetrate 
the  mortar  between  pieces  of  aggregate  forming  a  continuous  pattern 
appropriate  to  the  stress  condition.  Although  Hsu,  et  al .  (1978) 

studied  concrete  in  compression,  the  same  behavior  may  be  expected 
in  tension  but  developing  more  rapidly  than  in  compression.  Consider 
a  block  of  concrete  under  tension  in  direction  i  with  a  pattern  of 
internal  cracks  as  shown  in  Fig.  2.4.  Let  the  block  be  subjected 
to  tension  or  compression  in  any  direction  perpendicular  to 
direction  i.  The  contraction  or  expansion  in  direction  i  due  to 
the  applied  stresses  in  a  perpendicular  direction  will  be  absorbed, 
in  part,  by  filling  or  increasing  the  volume  of  cracks.  In  other 
words  the  total  Poisson's  ratio  effect  on  the  block  will  be  diminished. 
Once  the  material  has  reached  its  ultimate  strength  in  tension  in 
direction  i,  these  cracks  may  become  wide  enough  to  inhibit  inter¬ 
action  between  direction  i  and  the  normal  plane,  except  for  the  shear 
provided  by  aggregate  interlock,  etc.  Therefore,  it  has  been 
assumed  that  Poisson's  ratio  in  tension  is  constant  up  to  the 
appearance  of  mortar  cracks.  It  is  further  assumed  that  a  definite 
pattern  of  mortar  cracks  in  tension  starts  to  appear  at  50%  of  the 
ultimate  tensile  strength.  At  this  point  it  is  assumed  that  Poisson's 
ratio  starts  to  degrade  becoming  zero  at  ultimate  tensile  strength. 


■- 
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This  is  expressed  as  follows,  (see  Figures  2.5a  and  2.5b) 


V. 

i 


V 


o 


0  <  £  ./£  .  <  0.5  (2.34a) 

ui  ci  — 


V. 


l 


2V  (l-e  ./e  .) 

o  U1  Cl 


0.5  <£./£.  <  1.0  (2.34b) 

ui  ci  — 


V. 


0.0 


1.0  <  £  ./£  . 

UI  Cl 


(2.34c) 


l 


At  this  point  it  should  be  noted  that  the  variable  (J)  appearing 
in  Eq.  2.13  should  always  be  non  negative.  A  limiting  value  of  0.5 


has  therefore  been  placed  on  the  values  of  in  compression.  This 


value  corresponds  to  a  limit  of  zero  incremental  volume  change. 
Kotsovos  and  Newman  (1977)  have  noted  that  the  point  at  which  this 
limit  is  reached  corresponds  to  the  onset  of  unstable  micro-crack 
propagation.  This  is  the  process  that  causes  the  dilatation  phenomena 
observed  in  concrete  upon  approaching  the  ultimate  strength  under 
uniaxial  and  biaxial  compression.  Therefore,  it  can  be  argued  that 
the  dilatation  has  no  meaningful  constitutive  significance.  More¬ 
over,  no  dilatation  is  observed  under  triaxial  compression  (Schickert 
and  Winkler,  1977) .  Rather,  the  material  appears  to  flow  which 
supports  a  limit  of  0.5  on  Poisson's  ratio  in  compression. 

2.3.7  The  Ultimate  Surfaces 

The  evaluation  of  the  hypoelastic  moduli  described  in 

Section  2.3.5  and  2.3.6  requires  the  specification  of  the  parameters 

appearing  on  the  right  hand  side  of  Eqs .  2.21  to  2.34,  i.e.  O  •  and 

£  .  Since  these  parameters  vary  with  the  changing  stress  config- 

ci 

uration  as  well  as  the  loading  history,  the  evaluation  can  be  done 
by  specifying  a  surface  in  stress  space  to  define  the  three  values 


. 


23 


of  Q  .  and  a  surface  in  the  equivalent  uniaxial  strain  space  to 


cx 


define  the  three  values  of  £c^  that  correspond  to  the  cr^'s. 

A  surface  in  stress  space  that  defines  the  ultimate  strengths 

Q  .  for  any  ratio  of  stresses  is  called  a  "failure  surface".  This 
ci 

name  is  somewhat  misleading  in  case  of  a  material  with  strain 
softening  behavior.  It  is  proposed  to  call  such  a  surface  an 
"ultimate  strength  surface"  (Elwi  and  Murray,  1979) .  The  five 
parameter  surface  proposed  by  Wiliam  and  Warnke  (1975)  and  illu¬ 
strated  in  Fig.  2.6a  is  used  in  this  work  to  define  the  ultimate 
strength  surface  as  well  as  the  corresponding  equivalent  uniaxial 
strain  surface.  The  characteristics  of  this  surface  are  reviewed 
briefly  in  the  following. 


Let  the  mean  (average)  normal  stress  be  defined  as 


a 


(2.35) 


a 


and  let  the  mean  (average)  shear  stress  be  defined  as 


T 


a 


(2.36) 


where  s. .  is  the  deviatoric  stress  tensor  which  is  written  as 
13 


s 


ij 


(2.37) 


Dividing  by  the  uniaxial  compressive  strength  f  ,  the  variables 

Q  and  T  are  nondimensionalized  such  that 
a  a 


a 


(2.38a) 


a 


T 


a 


(2.38b) 
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The  surface  under  consideration  is  assumed  to  intersect  the 
deviatoric  plane  (a  =  constant)  in  three  symmetric  elliptical 

3. 

segments  forming  a  closed  convex  and  continuous  curve.  This  trace 
is  shown  in  Fig.  2.6b  on  a  triaxial  plot,  with  the  axes  representing 
the  main  deviatoric  stresses  0^ .  The  elliptic  trace  may  be  mathe¬ 
matically  described  as 

T  =  r (6,  0  )  (2.39) 

a  a 

where  0  is  called  the  angle  of  similarity  and  is  expressed  in  terms 
of  the  principal  stresses  as  (Wiliam  and  Warnke ,  1975) 

cos0  =  (ai  +  o2  -  2a3)/{/2[  (ai  -  a2) 2  +  (oz  -  03) 2  +  (a3  -  ai) 2]  ^2} 

(2.40) 

For  CJi  >_  O  2  ^  3  /  then  0  <_  0  <_  60° ,  as  may  be  seen  in  Fig.  2.6b.  The 

function  r  in  Eq.  2.39  is  defined  as  (Wiliam  and  Warnke,  1975) 

2  2  2  V2 

2rz^i  cos0  +  rz  (2ri  -  rz)  (4ri2cos  0+  5^  -  4rir2) 

r  =  -  (2.41) 

4ri2cos20  +  (r2~2ri)2 


in  which 


ri  2 


(2.42) 


The  variables  r2  and  ri  are  respectively  the  maximum  (0  =  60°)  and 
minimum  (0  =  0°)  radii  of  the  deviatoric  trace  of  the  surface.  These 
variables  are  assumed  to  be  parabolic  functions  of  the  hydrostatic 
stress  and  are  expressed  as  (Wiliam  and  Warnke,  1975) 


ri  =  a  +  a,  O  +  a9  O 
1  o  1  a  2  a 


(2.43a) 


r2  =  b  +bi  O  +b2  O 
2  o  a  a 


2 


(2.43b) 
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Plotted  on  the  so  called  "rendulic  plane",  the  functions  rl  and  r2 

are  illustrated  in  Fig.  2.6c.  The  values  of  the  coefficients  a  to 

o 

b2  in  Eqs .  2.43  are  chosen  such  that  the  variables  r:  and  r2  pass 
through  a  set  of  control  points,  as  illustrated  in  Fig.  2.6c. 

Appendix  A  describes  how  these  coefficients  are  evaluated.  The 
surface  proposed  by  Wiliam  and  Wanke  (1975)  is  now  completely  defined 
and  serves  to  evaluate  the  three  O  . 's  required  in  Sections  2.3.4 


and  2.3.5. 


However,  Eqs.  2.21  to  2.34  also  require  the  evaluation  of 

the  equivalent  uniaxial  strains  £  .  associated  with  <J  . .  For  this 

ci  ci 

purpose  it  is  postulated  that  there  is  a  surface  in  equivalent  uni¬ 
axial  strain  space  which  has  the  same  form  as  the  stress  surface 
described  above.  Therefore,  the  strain  quantities  and  y^  are 
defined  to  correspond  with  0 ^  and  respectively  by  replacing  , 
a2 ,  a 3  in  Eqs.  2.35  to  2.43  by  the  principal  components  of  the 
equivalent  uniaxial  strain  tensor.  The  analogous  nondimensionalization 
of  Eqs.  2.38  is  carried  out  with  £  ,  the  strain  corresponding  to 

uniaxial  compressive  strength,  replacing  f  .  Analogous  quantities 
for  the  strain  control  points  follow  directly  and  Eqs.  2.40  to 
2.43  together  with  the  equations  of  Appendix  A  can  be  used. 


2.3.8  Failure  Modes 

The  usage  of  the  ultimate  surfaces  as  implemented  by  Darwin 
and  Pecknold  assumes  that  instantaneously,  the  loading  is  pro¬ 
portional.  Thus,  joining  the  origin  and  the  current  stress  point 
with  a  straight  line  and  extending  the  straight  line  to  intersect 
the  surface  the  required  ultimate  point  can  be  obtained.  This 
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procedure  has  been  found  to  be  quite  adequate  in  triaxial  and 
biaxial  compression,  (Elwi  and  Murray,  1979)  .  However,  in  other 
situations  it  becomes  inadequate  since  it  is  conceivable  that  a 
point  under  tension  may  crack  normal  to  the  direction  of  loading, 
but  can  still  support  increasing  stress  in  the  other  two  directions 
in  tension  or  compression,  (Chitnuyanondh ,  et  al.,  1979).  Therefore, 
the  failure  mode  must  be  taken  into  consideration. 

For  the  purpose  of  imposing  failure  modes  on  the  evaluation 
of  the  parameters,  (J  ^  and  £c^ ,  a  number  of  conditions  are  stipulated 
as  follows : 

(i)  Compression-Compression-Compression : 

In  this  case  the  parameters  are  obtained  from  the  ultimate 
surfaces  described  in  Section  2.3.7. 

(ii)  Compression-Compression -Tension : 

In  this  case  the  tensile  components  are  obtained  from  the 

ultimate  surfaces.  The  compression  components  if  greater  than  f 

in  case  of  stresses  or  greater  than  £  in  case  of  equivalent 

uniaxial  strains  may  be  taken  equal  to  f  or  £  as  the  case  may 

cu  cu 

be.  The  term  "greater  than"  in  the  last  sentence  is  meant  in  the 
algebraic  sense. 

(iii)  Tension-Tension-Compression : 

In  this  case,  the  larger  tensile  components  are  obtained 

from  the  ultimate  surfaces.  The  smaller  tensile  components  may  be 

taken  as  f^  or  £^  as  the  case  may  be.  The  compressive  components 

may  be  taken  as  f  or  £  as  the  case  may  be. 

cu  cu 

(iv)  Tension-Tension-Tension: 

In  this  case  all  components  may  be  taken  as  f  or  £  as 
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the  case  may  be . 

The  stress  paths  that  may  produce  conditions  (ii)  to  (iv) 
are  illustrated  in  Figures  2.7a  to  2.7c. 

2.3.9  Verification  of  the  Constitutive  Relation 

In  order  for  the  constitutive  relation  proposed  in  this 
chapter  to  have  validity,  it  should  be  able  to  predict  strains 
associated  with  test  results  in  the  literature.  The  available 
experimental  studies  are  associated  with  a  predefined  stress  path. 

In  this  type  of  application  stress  increments  are  applied  and  it  is 
required  to  find  the  associated  strain  increments.  A  flow  chart  for 
this  type  of  application,  which  requires  iteration,  is  shown  in 
Fig.  2.16. 

A  comparison  of  predicted  strains  with  selected  sets  of 
measurements  by  Kupfer,  Hilsdorf  and  Rusch  (1969) ,  hereinafter 
refered  to  as  the  KHR  data,  in  a  series  of  biaxial  tests  is  shown 
in  Figs.  2.8  to  2.11,  inclusive.  The  parameters  that  have  been  used 
to  model  this  material  are  given  in  Table  2.1.  Fig.  2.8  shows  the 
uniaxial  compression  comparison,  and  Fig.  2.9  shows  the  biaxial 
compression  comparison.  Since  the  peak  stresses  and  corresponding 
strains  are  control  points  on  the  associated  ultimate  surfaces , 
exact  correlation  is  obtained  at  these  points.  Fig.  2.10  shows 
the  comparison  for  biaxial  compression  in  the  ratio  of  -l:-52. 

Since  the  chosen  ultimate  strength  surface  does  not  predict  the 
same  strength  as  the  KHR  failure  surface  there  is  a  discrepancy 
of  approximately  8%  in  this  instance.  A  comparison  of  typical 
response  in  tension-compression  is  shown  in  Fig.  2.11. 


. 
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The  above  comparisons  are  those  of  the  three  dimensional 
theory  with  plane  stress  tests.  A  comparison  of  the  theory  with 
some  of  the  three  dimensional  tests  of  Schickert  and  Winkler  (1977) , 
hereinafter  refered  to  as  the  SW  data,  has  also  been  carried  out. 
First,  the  parameters  of  Eqs .  2.43  have  been  evaluated  for  the  control 
points  and  are  shown  in  Table  2.2.  The  correspondence  of  these 
equations  on  the  rendulic  plane  with  the  SW  data  is  shown  in  Fig. 

2.12.  The  values  of  rx  and  r2  are  seen  to  agree  suitably  with 
this  data  set.  The  SW  tests  were  carried  out  by  loading  along  a 
hydrostatic  stress  path  with  deviatoric  stress  equal  to  zero,  and 
then  incrementing  the  three  principal  stresses  in  such  a  way  that 
the  hydrostatic  stress  remained  constant  while  the  deviatoric 
stress  varied.  Fig.  2.13  illustrates  results  for  deviatoric  stress 
path  along  line  0  -  A  of  Fig.  2.6b  (0  =  60°).  Fig.  2.14  illustrates 
results  for  a  deviatoric  stress  path  along  line  0  -  C  of  Fig.  2.6b 
(0  =  30°).  Fig.  2.15  illustrates  results  for  a  deviatoric  stress 
path  along  line  0-B  of  Fig.  2.6b  (0  =  0°)  .  It  can  be  seen  that 
the  model  reasonably  simulates  the  SW  data  along  each  of  these  non¬ 


trivial  stress  paths. 
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Ultimate  Strength 
Surface  Parameters 


Equivalent  Uniaxial 
Strain  Surface  Parameters 


Elastic  Moduli 


f  * 

-4650.0 

£ 

-0.00215 

E.  * 

4.7  *  10 

cu 

cu 

0  1 

a 

1.15 

a 

1.7512 

E  * 

4.7  *  10 

C 

c 

0  2 

a 

0.091 

a 

0.046 

E  * 

0 

t 

t 

0  3 

5, 

13.50 

Si 

50.0 

G  * 

0  1  2 

2.1  *  10 

Pi 

0.0 

Pi 

0.0 

Voi 

0.195 

4.0 

s2 

4.3 

vo2 

0 

P2 

0.0 

p2 

0.0 

V0  3 

0 

6 

6 

6 


Table  2.1  KHR  Material  Model  Parameters 
(*  in  psi) 


Ultimate  Strength 
Surface  Parameters 

Equivalent  Uniaxial 

Strain  Surface  Parameters 

Initial  Elastic 
Moduli 

f  * 

-4435.0 

£ 

-0.00283 

En  ,  * 

3.8  * 106 

cu 

cu 

0  1 

a 

1.21 

a 

1.3 

E  * 

3.8  *  106 

C 

c 

0  2 

at 

0.1 

at 

0.0461 

E03* 

3.8  *  106 

5, 

13.50 

Si 

22.5 

goi2* 

1.6  *  10G 

Pi 

0.0 

Pi 

0.0 

V0  1 

0.18 

$1 

3.75 

£ 2 

27.5 

V0  2 

0.18 

P2 

0.0 

P2 

0.0 

V0  3 

0.18 

Table  2.2  Schickert-Winkler  Material  Mode 
Parameters 
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Fig.  2.1  Compressive  Stress  Equivalent  Uniaxial 
Strain  Relation 


Fig.  2.2  Tensile  Stress  Equivalent  Uniaxial 
Strain  Relation 


Internal  Mortar  Cracks 
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Fig.  2.4  Mortar  Cracks 


(a)  Compression 


Fig.  2.5  Proposed  Poisson's  Ratio-Equivalent  Uniaxial 
Strain  Relation 
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(b)  A  Deviatoric  View  of  Argyris  Surface 
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(c)  The  Rendulic  View  of  Argyris  Surface 


Fig.  2.6  The  Ultimate  Strength  Surface 
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(a)  Compression-Compression-Tension 


t 


(c)  Tension-Tension-Tension 


Fig.  2.7  Failure  Conditions 
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STRRIN 


Fig.  2.8  Uniaxial  Compression  Comparison  with 
KHR  Data 


STRRIN 


Fig.  2.9  Biaxial  Compression  ( — 1 ;  —  1 )  Comparison 
with  KHR  Data 
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STRAIN 

Fig.  2.10  Biaxial  Compression  (— 1:— .52)  Comparison 
with  KHR  Data 


CM 


STRAIN 

Fig.  2.11  Compression  Tension  (-1:.204)  Comparison 
with  KHR  Data 


The  Compressive  Branch  (r, 
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Fig.  2.12  Rendulic  Plane  Comparison  with  Schickert  Winkler  Data 
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STRRIN 


.  2.13  Triaxial  Compression  Comparison  (-1:.5:.5 
stress  changes)  with  Schickert  and  Winkler 
Data 


STRAIN 


Fig.  2.14  Triaxial  Compression  Comparison  (-1:0:1 
stress  changes)  with  Schickert  and 
Winkler  Data 


SIG1/FCU 


STRRIN 

Fig.  2.15  Triaxial  Compression  Comparison  (.5:  .5:— 1 

stress  changes)  with  Schickert  Winkler  Data 


Fig.  2.16  Flow  Chart  for  a  Specified 


Stress  Path 


CHAPTER  THREE 


FINITE  ELEMENT  MODEL 


3. 1  Introduction 

As  described  in  Chapter  One,  it  is  the  purpose  of  this  study 
to  develop  a  sophisticated  tool  for  the  analysis  of  reinforced  con¬ 
crete  structures,  for  axisymmetric  and  planar  problems,  under  static 
short  term  loads  and  small  displacement  fields.  The  rectangular 
serendipity  family  of  finite  elements  (Zienkiewicz ,  1971)  has  been 
chosen  as  the  basis  of  program  FE PARCS 5 .  (FEPARCS  is  an  acronym  for 

'Finite  Element  Program  for  the  Analysis  of  Reinforced  Concrete 
Structures'.)  These  elements  can  be  converted  easily  from  plane 
stress  to  axisymmetric  state  or  vice  versa  without  changing  the  frame 
of  reference  or  the  description  of  the  strain  field  (Zienkiewicz, 
1971) .  The  main  difference  between  axisymmetric  and  planar  problems 
is  in  the  description  of  material  properties  and  the  addition  of  cir¬ 
cumferential  strains  (in  case  of  axisymmetric  formulation) . 

Three  types  of  elements  were  developed;  the  solid  element; 
a  meridional  reinforcing  element  and  a  circumferential  reinforcing 
element.  The  meridional  reinforcing  element  represents  the  rein¬ 
forcing  bars  or  prestressing  tendons  acting  in  z  -  r  plane.  The 
z-axis  is  the  vertical  axis  of  symmetry  and  the  r-axis  is  horizontal 
radial  axis.  The  circumferential  reinforcing  element  represents 
the  circumferential  reinforcing  bars  or  prestressing  tendons.  These 
elements  are  empty  solid  elements  except  for  layers  of  reinforcing 
or  prestressing.  Strain  compatibility  between  the  steel  and  the 
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concrete  is  preserved  by  using  the  same  shape  functions  and  order  of 
integration  for  all  elements,  (Buyucozturk ,  1977). 

A  variational  approach  was  used  to  derive  the  element 
stiffness  matrices  as  well  as  different  types  of  element  loads  such 
as  gravity,  thermal,  and  psuedo  loads  and  surface  tractions.  In 
this  chapter  an  incremental  variational  principle  is  reviewed, 
the  element  stiffness  matrices  are  formulated,  the  loads  resulting 
from  body  forces  and  surface  tractions  are  derived  and  the  boundary 
conditions  are  described. 

3 . 2  Fundamental  Variational  Formulation 

Before  attempting  to  formulate  the  various  element  property 

matrices  and  vectors,  it  is  necessary  to  establish  the  functional 

relationships  governing  the  problem.  The  variational  approach 

provides  a  powerful  formulation  technique  for  finite  element  theory. 

In  this  section  the  necessary  variational  principles  are  presented 

in  incremental  form  based  on  small  displacement  fields. 

Let  the  body  shown  in  Fig.  3.1  be  divided  into  "K" 

number  of  elements.  Let  V  be  the  volume  of  an  element  and  let  S 

e  qe 

and  S  be  those  portions  of  the  element  surface  on  which  the  dis- 
Oe 

placement  and  stresses  are  prescribed  respectively.  The  superscript 
(  ) 0  indicates  initial  quantities  at  the  beginning  of  a  load  step, 
and  the  prefix  A(  )  denotes  incremental  quantities.  Prescribed 
quantities  will  be  denoted  by  (  ) .  The  incremental  field  equations 


can  be  written  as 


. 
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a? .  .  +  Ag. .  .  +  f. 0  +  Af.  =  0 

13/3  13  g  1  1 


(3.1a) 


Ae.  . 

id 


4  (q?  .  +  q°  )  +4  (iq.  .  +  Aq .  .)-£?. 
2  I'D  3  1  2  1 ' D  D'1  ID 


(3.1b) 


aaij  ■  cijkJiAe^ 


(3.1c) 


Eqs.  3.1a  to  3.1c  are,  respectively,  the  equilibrium  equation,  the 
strain  displacement  relation  and  the  instantaneous  stress  increment 
strain  increment  relation.  The  symbols  G,  e,  F,  q  and  C  denote 
respectively  the  stress  tensor,  the  strain  tensor,  the  body  force 
per  unit  volume  vector,  the  displacement  field  and  the  constitutive 
tensor  (Pian,  1976) . 

The  mechanical  boundary  conditions  on  are  defined  as 

(Pian,  1976) 


T.  0 

1 

=  G° .  n . 

ID  D 

(3.2a) 

At. 

1 

=  Aa . .  n . 

ID  D 

(3.2b) 

in  which  T  denotes  surface  tractions  and  n  denotes  a  unit  vector 

normal  to  the  surface.  The  displacement  boundary  conditions  on 

S  are  written  as  (Pian,  1976) 
qe 

q. 0  +  Aq.  =  q. 0  +  Aq.  (3.3) 

1  1  1  1 

Writing  the  first  variation  of  Eqs.  3.1b  and  3.3  and  noting 
that  initial  and  prescribed  quantities  do  not  vary,  the  following 
conditions  are  obtained. 

=  (6  Aq.  .  +  6  Aq .  .)  (on  V  ) 

2  1 1  j  3  / 1  e 


6  Ae.  . 
ID 


(3.4a) 
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6  Aq .  =  0 

l 


(  S  ) 


qe 


(3.4b) 


Let  the  structure  be  forced  through  a  virtual  displacement 
field  5  Aq  subject  to  the  condition  stated  in  Eq.  3.4b.  The  virtual 
work  6u  of  the  internal  forces  and  the  virtual  work  6w  of  the  external 
forces  are  expressed  as 


5u 


£  { 
K 


(a°  .  +  Aa.  .)  SAe.  .  dv} 

id  id  id 


(3.5a) 


Sw 


£  { 
K 


+  A  F.  ) 

1 


dV  + 


(T.°  +  A  T.  ) 

l  l 


6Aq_^  ds} 


(3.5b) 


Substituting  Eqs .  3.2  into  the  second  term  of  Eq.  3.5b,  applying 
the  divergence  theorm  to  the  same  term  and  using  Eqs.  3.1a,  3.4a 
and  3.4b  it  can  be  proven  that  the  virtual  work  of  internal  forces 
equals  the  virtual  work  of  external  forces  subject  to  the  conditions 
stated  in  Eqs.  3.4  which  is  called  the  "principle  of  virtual  work" 
and  can  be  written  as 


£  { 
K 


[ (a0 .  +  Aa. .)  6Ae.  . 

id  id  id 

V 

e 


(F?  +  AF. )  6Aq . ]  dV 

ill 


(T  0  +  At. )  5Aq .  ds}  =  0  (3.6) 

l  l  1  1 

e 


If  the  initial  state  is  in  equilibrium,  then  the  virtual 
work  associated  with  it  must  vanish.  In  this  case  Eq.  3.6  reduces 
to  the  form 


^  { 
K 


•V 


[Aa.  .  6Ae..-Af,  6Aq ,  ]  dv  - 

ID  ID  i  i 


At.  6Aq .  ds}  = 

i  i  J 


0 


(3.7) 
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which  is  the  incremental  form  of  the  principle  of  virtual  work 
(Horrigmoe  and  Bregan,  1976)  . 

The  potential  energy  7T^  associated  with  an  increment  of 
displacement  Aq  is  defined  as  (Pian,  1976) 


7T  (Aq)  = 
P 


Cijk£Aek£ 


(F?  +  AF. )Aq. ]  dV 

l  l  i 


(T.°  +  At.0  )  Aq.  ds}  (3.8) 

l  l  l  J 

e 


Carrying  out  the  first  variation  on  Eq.  3.8  the  right  hand  side 
becomes  identical  to  the  left  hand  side  of  Eq.  3.6.  This  means 
that  the  potential  energy  is  stationary  subject  to  the  conditions 
stated  in  Eqs.  3.4.  It  is  understood  here  that  the  material  is 
linear  elastic  in  the  incremental  region.  In  addition  it  can  be 
proven  that  if  the  strain  energy  denoted  by  (y  Ae^_.  Ae^) 

positive  definite,  then  the  stationary  value  of  the  potential 
energy  is  a  local  minimum  (Fung,  1965) .  These  two  principles 
can  be  mathematically  stated  as 


6tt  =  0  (3.9a) 

P 

7T  (Aq  +  6Aq)  -  TT  (Aq)  >  0  (3.9b) 

P  P  - 

Having  presented  the  incremental  variational  principles  it 
becomes  possible  to  formulate  the  set  of  equations  necessary  to  solve 


the  problem.  Eq.  3.8  is  written  in  matrix  form  as 


. 
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7T 

P 


[<  a  >°{  Ae  }  -  <  F  >°{  Aq  }]dv 


* 

<  T°  >{  Aq  }  ds 
e 


+ 


[y  <  Ae  >  [  C  ]  {  Ae  }  -  <  Af  >{  Aq  }  ]  dv 
v  2 

e 


* 

<  At  >{  Aq  }  ds}  (3.10) 

•'s 

e 


Using  interpolating  functions,  the  displacement  and  strain  fields  can 
be  described  in  terms  of  the  nodal  displacements  as 


{Aq}  = 

[  N  ]{  Aq  } 

(3.11a) 

{Ae}  = 

[  B  ]  {  Aq  } 

(3.11b) 

where  [N] 
operators 
Eqs .  3.11 


is  the  matrix  of  interpolating  functions,  [B]  is  a  differential 
matrix  and  {Aq}  are  the  nodal  displacements.  Substituting 
into  Eq.  3.10,  the  latter  assumes  the  form 


TT 

P 


I  { 

K 


[<  G°  >[B]  {  Aq  }  -  <  F°  >[N]  {  Aq  }  ]  dV  - 


V 


<  T°  >  [N]  {  Aq  }  dS 


+ 


[-|  <  Aq  >[B]T[C]  [B]{  Aq  }  -  <  AF  >[N]{  Aq  }]  dV 
V 

e 


<  At  >[N]{  Aq  }  ds}  (3.12) 

J  S 

e 

The  summation  over  the  elements  making  up  the  continuum  under 
consideration  can  be  carried  out  over  the  individual  terms  of 


Eq.  3.12  yielding 


■ 
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77  =  j  Z  <  Aq  >  [  K  ]  {  Aq  }  -  E  <  Aq  >{  AF  } 

M  ^  0  0 


Z  <  Aq  >{  At  1  +  <  Aq  >  {  Aq  } 

K  e  e 


(3.13) 


where 


tKe] 


V 


[  B  ]  [  C  ]  [  B  ]  dV 


(3.14a) 


{  Af  } 
e 


V 


T  — 

[  n  ]  {Af  }  dv 


(3.14b) 


{  At  } 

e 


[  N  ] T  {  AT  }  dS 


(3.14c) 


{  Aq  } 

e 


[  B  ]T  {  CT°  }  -  [  N  ]'  {  F  0  }]  dV- 


V 


T;„-o 


[  N  ]T  {t0}  dS 


(3 . 14d) 


The  terms  [K  ]  ,  <  Af  >  and  <  At  >  represent,  respectively,  the 
e  e  e 

element  stiffness  matrix,  the  increment  of  body  forces  and  the 

increment  of  surface  tractions.  The  term  <  Ao  >  should  vanish  if 

e 

the  initial  state  is  in  equilibrium  and  may  represent  the  unbalanced 
forces  otherwise.  The  principle  of  Eq.  3.9a  can  be  applied  to 
Eq.  3.13  subject  to  the  conditions  stated  in  Eqs .  3.4,  recognizing 
the  arbitrary  nature  of  <6Aq>  ,  yielding 


E  [  k  ]  {  Aq  >  -  E  {  Af  }  -  Z  {  At  }  +  Z  {  AQe> 

K  6  K  6  K  S  K 


=  0 


(3.15) 


It  must  be  noted  that  the  summations  of  Eq.  3.15  are  carried 


out  in  the  direct  stiffness  sense  and  can  be  summarized  as 


[  K ]  {  Ar  }  -  (  Ar  }  +  {  AQ  } 


0 


(3.16) 
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where  [ K ]  is  the  structure  stiffness  matrix,  <  Ar  >  is  the  incre¬ 
ment  of  nodal  displacements,  <  Ar  >  is  the  increment  of  prescribed 
loads  and  <  A Q  >  is  the  unbalanced  load  at  the  end  of  the  previous 
load  step.  Eq.  3.16  together  with  the  condition  of  Eq.  3.4b  make 
up  a  nonsingular  system  of  equations  which  can  be  solved  for  the 
increment  of  displacement  <Ar>. 

3 . 3  Isoparametric  Element  Formulations 
3.3.1  Solid  Element  Formulation 

The  formulation  of  a  finite  element  depends  on  the  unique 
description  of  unknown  functions  within  the  element  in  terms  of 
parameters  associated  with  the  values  of  the  functions  at  the  nodes, 
(Zienkiewicz ,  et  al.,  1970).  These  parameters  are  based  on  what  are 
called  "shape  functions" .  The  shape  functions  depend  on  the  spatial 
description  of  the  element  and  must  satisfy  two  conditions: 

a)  Continuity  of  the  function  and  derivatives  up  to  the 
order  required  by  the  variational  formulation  which  is 
one  order  less  than  the  maximum  order  of  differentiation 
appearing  in  the  functional  (continuity  of  class  C°  is 
required  in  this  study) . 

b)  The  ability  to  represent  a  constant  field  of  the  unknown 
function  over  the  entire  volume  of  the  element. 

Condition  (a)  can  be  satisfied  if  the  number  of  nodes  on  a  boundary 
is  sufficient  to  determine  uniquely  the  variation  of  the  function 
along  the  boundary  and  if  only  the  interboundary  nodes  between  two 
adjacent  elements  influence  the  value  of  the  unknown  function  on 


. 
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this  boundary.  In  this  study  only  linear,  quadratic  and  cubic  rect¬ 
angular  isoparametric  elements  of  the  serendipity  family  are  used 
(see  Fig.  3.2).  The  shape  functions  for  these  elements  are  given 
in  Appendix  B. 

Since  the  elements  used  are  isoparametric  elements,  the  same 
shape  functions  can  be  used  to  interpolate  geometry  as  well  as  dis¬ 
placements  as  follows 

<  r  z  >  =  <<|>>[  _r  _z  ]  (3.17a) 

<uv>  =  <(J)>[_uv]  (3.17b) 

where  r  and  z  are  horizontal  radius  and  vertical  coordinate  respec¬ 
tively  (axisymmetric  formulation) ,  u  and  v  are  the  horizontal  and 
vertical  displacements  respectively,  the  symbol  (_)  denotes  nodal 
values  and  <  (f)  >  is  the  set  of  shape  functions. 

A  differential  element  of  volume  dV  representing  a 
radian  sector  in  any  axisymmetric  continuum  can  be  written  as 

dV  =  rdA 


where  dA  is  an  element  of  area  and  is  defined  as  the  magnitude  of 
the  cross  product  of  the  two  vectors  dbi  and  db2  making  up  the  sides 
of  the  element  of  area  as  shown  in  Fig.  3.3b.  Writing  the  cross 
product  in  terms  of  the  r  and  z  components,  Eq.  3.18  assumes  the  form 


dV 


dr  i 

dz  i 

dr2 

dz2 

(3.19) 


The  components  dr  and  dz  can  be  written  in  terms  of  the  normalized 
coordinates  E,  and  y  (see  Fig.  3.3a)  as 
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<  dr  dz  > 


[  J]  = 


=  < 

d£  dy  >  [  j  ] 

(3.20a) 

~9r 

9z 

3C 

3C 

(3.20b) 

9r 

9z 

9y 

9y 

where  [ J ]  is  the  Jacobian  matrix.  This  matrix  is  expressed  in 
terms  of  the  shape  function  derivatives  as 


[r  z]  (3.21) 

Since  dA  is  arbitrary  it  can  be  chosen  such  that  dbi  corresponds  to 
<  d£  o  >  and  db2  corresponds  to  <  o  dy  >  .  Thus  Eq.  3.19  can  be 
written  as 


[  J]  = 


<  -^  > 


3? 


<  4*  > 


dV  =  r  j  J  |  d^dy 


(3.22) 


where  the  symbol  (|  |)  denotes  a  determinant. 

The  purpose  of  this  section  is  to  evaluate  the  integral 
of  Eq.  3.14a.  Assuming  that  the  constitutive  matrix  [ C  ]  is  known 
and  that  the  element  of  volume  can  be  represented  by  the  right  hand 
side  of  Eq.  3.22,  it  remains  to  evaluate  the  matrix  [  B  ]  in  terms  of 
the  shape  functions  and  derivatives.  Let  the  strains  of  Eq.  3.11b 
be  expressed  for  axisymmetric  strain  fields  as  (Zienkiewicz ,  1971). 


— 

Ae 

r 

<  9(J)/9r  > 

0 

A  e 

z 

0 

<  9(f)/9r  > 

<  <j p/r  > 

0 

Ay 

•  rz 

<  9cJ)/9z  > 

<  9(|)/3r  > 

{Au} 


{A  v} 


(3.23) 


• ' 

1 
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The  derivatives  appearing  in  Eq.  3.23  may  be  evaluated  as 


<  9<j)/9r  > 

_  i 

<  3<f>/9£  > 

<  d§/dz  > 

=  [  J] 

<  3(J)/3y  > 

(3.24) 


Since  the  matrix  [  J  ]-1  has  polynomials  in  the  denominator 
it  becomes  impractical  to  try  to  integrate  Eq.  3.14a  in  closed  form. 
Therefore,  integrations  are  evaluated  numerically.  For  the  purposes 
of  this  study,  Gaussian  integration  was  selected  as  a  viable  scheme. 
Thus  Eq.  3.14a  may  now  be  written  as 


[K 


n  m 

I  I  [b(£. ,y .) ]t[c] [b(£. ,yi) ]wnwm  |j|r 

i=l  j=i  1  3  1  1  J 


(3.25) 


where  n  and  m  are  the  order  of  Gaussian  integration  in  directions 
E,  and  y  respectively,  and  w^  and  w™  are  the  weights  associated 
with  each  point  in  directions  E,  and  y  respectively. 


3.3.2  Formulation  of  the  Meridional  Reinforcing  Isoparametric  Element 
In  order  to  model  the  meridional  reinforcing  and  prestressing 
layers  within  a  solid  element,  a  compatible  element  can  be  described 
such  that  the  integration  is  carried  out  only  along  the  required 
layers.  An  element  which  can  accommodate  layers  along  the  y  or  E, 
directions  is  formulated  in  this  section,  see  Fig.  3.4a. 

Let  the  potential  energy  of  the  element  be  written  as 


IT  (Ae) 
P 


(a°Ae  +  y  AaAe)  dv 
•  v 


(3.26) 


e 

where  the  strain  and  the  stress  are  defined  only  along  the  layer 
under  consideration.  The  element  of  volume  dV  is  written  in  terms 


‘ 
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of  the  element  of  length  along  the  layer  d£,  the  area  of  the  layer 
per  unit  width  A  and  the  radius  (axi symmetric  formulation)  as 

dV  =  r  Ad£  (3.27) 


The  strain  increment  is  expressed  in  terms  of  the  global  strain  and 
the  orientation  of  the  layer  as 

Ae  =  Ae  cos20  +  Ae  sin20  +  Y  sin0cos0  (3.28) 

r  z  rz 


where , 


COS0 


dr 

d£ 


sin0 


dz 

d£ 


(3.29a) 


(3.29b) 


Using  the  definition  of  the  strain  increment  stated  in  Eq.  3.1b 
neglecting  initial  strains  and  using  Eqs .  3.29,  Eq.  3.28  may  be 
written  as 


Ae 


3Au  dr  dr  3Av  dz  dz  3Au  dr  dz  3Av  dr  dz 
3r  d£  d£  3z  d£  d£  3z  d£  d£  3r  d£  d£ 


Eq.  3.30  can  be  rearranged  such  that 


dAu  dr  |  dAv  dz 
d£  d£  d£  d£ 


(3.31) 


Let  the  layer  under  consideration  be  in  the  y  direction. 
The  element  of  length  d£  and  the  direction  cosines  are  expressed 
in  terms  of  the  normalized  coordinates  as 


d£  =  J'  dy 


-  —  /  J' 

“  3y  7 


(3.32a) 


COS0 


(3.32b) 
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sin0  =  /  J' 

9y 


(3.32c) 


where  J  is  a  transformation  factor  which  can  be  derived  using  Eqs .  3.20 
and  is  written  as 

1/2 


J'  = 


(9r/8y) 2  +  (3z/3y) 2 
Eq.  3.31  can  now  be  stated  as 


.  _  dAu  dy  dr_  dy_  dAv  dy  dz  dy 

dy  d£  dy  d£  +  dy  d£  dy  d£ 


(3.33) 


(3.34) 


Using  the  definitions  of  Eqs.  3.32,  Eq.  3.34  can  be  simplified  as 


Ae  = 


dAu  dAv 

cost)  +  — ; — sint) 


Idy 


dy 


/  J' 


(3.35) 


The  increment  of  stress  appearing  in  Eq.  3.26  can  be  stated 
in  terms  of  the  increment  of  strain  Ae  and  the  current  elastic 
modulus  E  as 


Aa  =  e  Ae 


(3.36) 


Substituting  for  A<J  and  dV  in  Eq.  3.26  using  Eqs.  3.27  and  3.36 
yields 


Up  = 


[a°  Ae  +  j  Ae  e  Ae]  rAdil 


(3.37) 


v 


The  second  term  of  Eq.  3.37  yields  the  element  stiffness  matrix  as 


COS0 

r~ 

o 

A 

^1 

z  { 

I-1 

E  <  cos0  sin0  > 

rS  1 

-e- 

o 

_ 1 

sin0 

1 

O 

A 

n 

CL 

Ar  itk 

^wi} 


(3.38) 
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where  wm  is  the  Gaussian  weight  and  m  is  a  suffix  denoting  the  order 
of  Gaussian  integration.  When  the  layer  is  along  the  E,  direction,  E, 
is  substituted  for  y  in  all  the  pertinent  equations  of  this  section. 

3.3.3  Formulation  of  A  Circumferential  Reinforcing 
Isoparametric  Element 

Hoop  reinforcing  and  prestressing  layers  can  be  modelled  in 
the  same  manner  as  the  meridional  layers  described  in  Section  3.3.2. 
Again,  an  element  which  can  accomodate  circumferential  layers  parallel 
to  the  ]i  and  E,  coordinates  has  been  formulated.  In  this  case  the 
potential  energy  of  the  element  is  defined  as  stated  in  Eq.  3.26.  The 
element  of  volume  is  defined  as 

dV  =  rAd£  (3.39) 

where  A  is  the  area  of  the  layer  under  consideration  and  d£  is  an 
element  of  length  parallel  to  the  direction  of  the  layer.  The  strain 
however  is  defined  as 

Ae  =  Au/r  (3.40) 

Let  the  layer  under  consideration  be  in  the  y  direction.  The 
element  of  length  d£  is  defined  in  terms  of  the  normalized  coordinate 
y  as 

d£  =  J' dy  (3.41) 

where  J  is  the  transformation  factor  defined  by  Eq.  3.33.  For 
this  element,  the  stiffness  matrix  is  obtained  from  Eq.  3.37  by 
substituting  the  right  hand  side  of  Eq.  3.36  for  the  increment 
of  stress  and  Eq.  3.40  for  the  strain  increment  as 
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[Ke] 


=  z  { 


{(f)} 


E  <  <  d)  >  <  o  >  >  A  —  Wm} 

r  1 


(3.42) 


where  w™'s  are  the  Gaussian  weights  and  m  denotes  the  Gaussian  order 


of  integration. 


3 . 4  Body  Forces 

The  body  forces  that  arise  in  this  study  are  gravity  loads 
and  equilibrating  loads.  The  latter  are  equivalent  to  the  internal 
forces  in  the  structure.  Although  this  type  of  load  is  not  applied 
to  the  structure  its  evaluation  is  required  to  check  equilibrium  and 
to  evaluate  the  unbalance  between  the  applied  loads  and  the  state  of 
stress  in  the  structure.  In  this  section  gravity  loads  and  equili¬ 
brating  loads  as  well  as  thermal  loads  arising  from  possible  thermal 
fields  are  formulated  on  a  work  equivalent  basis. 


3.4.1  Gravity  Loads 

In  this  study  gravity  loads  are  assigned  to  the  negative  z 
direction  and  are  calculated  for  the  solid  element  alone.  Let  AF 
be  an  increment  of  gravity  load  per  unit  volume.  The  virtual  work 
done  by  this  load  in  going  through  the  virtual  displacement  6Aq,  which 
is  compatible  with  the  condition  stated  in  Eq.  3.4b,  is  the  second 
term  on  the  right  hand  side  of  Eq.  3.7  and  can  be  isolated  as 


6  (AF) 


y6Av  dV 


(3.43) 


where  y  is  the  specific  weight  of  the  solid  element  material.  In 
matrix  form  Eq.  3.43  can  be  written  as 


' 


-  <  <5Au  6Av  > 


dV 


(3.44) 


<  6^u  6Av  >  {  A  F  } 
—  —  e 


Y 

Ve 


m 


Substituting  for  the  element  of  volume  dV  in  Eq.  3.44  by  the  right 

hand  side  of  Eq.  3.22  and  carrying  out  the  integration  numerically, 

the  work  equivalent  gravity  loads  <  Af  >  are  expressed  as 

e 


{  A  F  ) 
e 


n  n 
-  Y  £  Z 
i=l  j=l 


T  n 

Jr  w. 

l 


m 
w . 
3 


(3.45) 


where  J,  r,  w,  n  and  m  have  been  defined  in  Section  3.3.1. 


3.4.2  Thermal  Loads 

Thermal  loads  arise  when  deformation  due  to  the  presence  of 
a  temperature  distribution  is  partially  or  fully  constrained  by  the 
presence  of  a  gradient,  external  indeterminate  restraints,  variations 
in  material  thermal  properties  or  internal  indeterminacy.  Let  the 
increment  of  strain  arising  in  a  body,  where  one  or  more  of  the  above 
described  situations  exist,  be  written  as 


A  eij  =  hjki  Aakji +  A£ij 


(3.46) 


where  A.  0  is  the  compliance  tensor  associated  with  C_.  „  defined 
13k*  ^  i3k£ 

in  Chapter  Two  and  Ae . is  the  increment  of  strain  occuring  due  to 
thermal  change  alone  and  is  defined  as 


A  =  a.  .  At  (3.47) 

13  13 

where  At  is  the  increment  of  temperature  and  a.  .  is  the  thermal 

13 

expansion  tensor.  It  must  be  noted  that  in  an  orthotropic  material 


shear  stresses  in  the  orthotropic  directions  do  not  occur  due  to 


■ 
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thermal  changes. 


In  the  principal  axes  of  orthotropy  a. .  is  defined  as 

iD 


a.  .  =  a. 

n  i 

(no  sum) 

(3.48a) 

a.  .  =  0 

13 

(i  *  j) 

(3.48b) 

where  a.  is  the  thermal  expansion  coefficient  in  direction  i. 
Inverting  and  rearranging  Eq.  3.46  the  increment  of  stress  due  to 
the  increment  of  temperature  can  be  written  as 


Ci  jk£ 


Ignoring  external  work, 


Eq.  3.7  is 


At) 

reduced  to  the  form 


(3.49) 


Z 

K 


<A  £k£  '  \l  At) 


c.  ..  0  6  A  e.  . 

lOkiC  13 


dV>  =  0 


V 


(3.50) 


The  first  term  of  equation  3.50  gives  rise  to  the  element  stiffness 
matrix.  The  second  term  gives  rise  to  the  thermal  loads.  Using 
Eq.  3.4a  and  Eq.  3.11b,  the  second  term  of  Eq.  3.50  is  stated  in 
matrix  form  as 


{AFt> 


z  { 

K  V 


[c  ] 


{a}'  At  dv} 


(3.51) 


e 

The  increment  of  temperature  At  is  expressed  in  terms  of  the  nodal 
values  and  shape  functions  as 


At  =  <  <f>  >  {  At  } 


(3.52) 


Using  Eq.  3.52  and  numerical  integration,  Eq.  3.51  can  be  written  as 


■ 
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(Af^} 


Z 

K 


n 

Z 


m 


m 


Z  [B(£.,y.)]  [c]  {a}<(})(£.,y  . )  >rJwn 


i=lj=l 


{At} 


(3.53) 


The  foregoing  formulation  has  been  for  a  solid  element. 
Similar  formulations  can  be  carried  out  for  the  meridional  and 
the  circumferential  reinforcing  elements  described  in  Sections  3.3.2 
and  3.3.3  resulting  respectively,  in 

(Af  }  =  Z 

K 


n 

r  4} 

0 

j  sin0 

■  n 

z 

i=V 

0 

{<!>}_ 

j  COS0 

■  <(t>>  OtAr  Ew. 

l 

{At}  (3.54a) 


{At}  (3.54b) 

where  the  appropriate  definitions  of  Sections  3.3.2  and  3.3.3  apply. 

The  summations  enclosed  in  square  brackets  in  Eqs.  3.53, 

3.54a  and  3.54b  are  carried  out  for  a  particular  solid  element  and 
all  the  meridional  and  circumferential  reinforcing  layers  within, 
summed  up  and  then  multiplied  once  by  the  element  nodal  temperatures. 
In  this  manner,  efficiency  of  calculations  is  increased.  Eqs.  3.53 
and  3.54  are  summarized  in  one  equation  as 

{At  }}  (3.55) 
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[Th]  +  [Th]  +  [Th] 
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<<(f)><o>>  EOtr  A  J  w ^ 


where  [Th]  ,  [Th]  and  [Th]  and  the  summations  enclosed  in  square 
s  mr  cr 


brackets  in  Eqs.  3.53,  3.54a  and  3.54b  respectively. 


58 


3.4.3  Equilibrating  Loads 

Equilibrating  loads  are  fictitious  loads  that  are  work  equi¬ 
valent  to  the  state  of  stress  in  a  structure.  The  estimation  of 
these  loads  is  necessary  to  check  equilibrium  and  to  calculate  the 
unbalanced  loads  in  case  equilibrium  is  not  satisfied.  The  equili¬ 
brating  loads  associated  with  an  element  constitute  the  first  integral 
on  the  right  hand  side  of  Eq.  3.14d.  Thus,  for  a  solid  element, 
using  Gaussian  integration  and  the  definition  of  the  element  of 
volume  in  Eq.  3.22,  the  equilibrating  loads  are  evaluated  as 

n  m 

{  Q  }  =  Z  Z  [B(£.,y.  )]  {0}°  r  J  wn  wm  (3.56) 

6  i=l  j=l  11  1  D 


The  equilibrating  loads  associated  with  a  meridional  reinforcing 
element  with  layers  running  in  the  y  and  £  directions  can  be 
evaluated  respectively  as 


m 

{  C>  }  =  Z 

rt*.u  HoJi 

sin0 
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i— i. 

ii 
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L_  — 

COS0 

(3.57a) 
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(4),^i}  {o} 

COS0 
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{o} 
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G  r  A  wi 


(3.57b) 


For  a  circumferential  reinforcing  element  the  equilibrating  loads 
are  written  as 


{  Q  } 
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=  Z 

i=l 
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{o} 
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A  J  w.  G 
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(3.58) 
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3 . 5  Surface  Tractions 

Program  FEPARCS5  accepts  surface  tractions  in  the  form  of 
nodal  pressure  intensities  which  may  be  normal  or  tangential  to  an 
element  surface.  An  element  surface  is  defined  as  a  group  of  adjacent 
nodes  forming  one  side  of  an  element.  The  normal  pressure  component 
is  positive  when  prescribed  along  the  right  hand  normal  to  an  element 
surface.  The  tangential  pressure  component  is  positive  when  it 
describes  90°  counterclockwise  angle  with  right  hand  normal  to  the 
element  surface,  see  Fig.  3.5a.  The  right  hand  normal  to  an  element 
surface  points  to  the  right  when  the  nodes  which  define  the  element 
surface  are  traversed  in  the  order  they  are  specified. 

The  third  integral  on  the  left  hand  side  of  Eq.  3.7  represents 
the  virtual  work  done  by  surface  tractions  when  acting  through  a 
virtual  displacement  field  which  satisfies  the  condition  stated 
by  Eq.  34b.  This  statement  is  mathematically  expressed  as 

* 

<  6Aq  >  {  At  }  =  E  <  5Aq  >  {  Ap  }  ds  (3.59) 

K  JsK 

where  <  Ap  >  is  the  pressure  field,  < 6Aq  >  is  the  virtual  displacement 
field  and  <  At  >  is  the  nodal  force  vector  equivalent  to  the  pressure 
field.  The  pressure  field  is  defined  in  terms  of  vertical  and 
horizontal  components,  or  alternatively  in  terms  of  the  normal  and 
tangential  components  as  shown  in  Fig.  3.5b  as  follows 

=  <  Ap  Ap  > 

r  z 


<  Ap  > 


(3.60) 


60 


sina  -cosa 

Apl 

{  Ap}  = 

n 

cosa  sina 

Apt 

(3.61) 


where  a  is  the  angle  which  the  right  hand  normal  to  the  surface 
describes  with  the  vertical  axis  as  shown  in  Fig.  3.5c,  and  is 
defined  mathematically  as 


sina  =  dz/d£ 

(3.62a) 

cosa  =  dr/d£ 

(3.62b) 

where  d£  is  an  element  of  length  along  the  surface.  In  terms  A  of 
a  normalized  coordinate,  e.g.  y ,  d£  can  be  written  as 

d£  =  J'dy  (3.63) 


where  J*  is  a  transformation  factor  defined  as 


/,9r  2  ,  ,3z  2 

=  ■'V  +  V 


(3.64) 


The  element  of  surface  dS  appearingin  Eq.  3.59  is  defined  for  a 
one  radian  sector  as 


dS  =  r  d  £  (3.65) 

Using  Eqs .  3.65,  3.63  and  3.61,  the  surface  tractions  are  expressed 
in  terms  of  shape  functions  and  normal  and  tangential  components  as 


{  At  }  =  Z 

,  1 

(4)} 

{o}“ 

sina  -cosa 

^>n 

K 

i 

_{o} 

(4>}_ 

cosa  sina 

rJ'  dp  (3.66) 


Using  numerical  integration  and  writing  the  pressure  components  in 


terms  of  nodal  pressure  intensities,  Eqs.  3.66  can  be  written  as 
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{  At  }  =  £  Z 
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o  { 4)  } 
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<  (f)  >  o 
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{  Ap  } 


-n 


(3.67) 


If  the  conditions  of  compatibility  between  elements  described 
in  Section  3.3.1  are  satisfied,  the  same  shape  functions  defined  in 
Appendix  B  can  be  used. 


3 . 6  Boundary  Conditions 

Boundary  conditions  in  program  FEPARCS5  are  applied  in  the 
form  of  linear  springs  of  very  high  stiffness  compared  to  the 
stiffness  of  the  structure  under  consideration.  These  springs  can 
have  any  orientation  in  the  plane  r  -  z .  Let  the  virtual  work 
associated  with  a  boundary  element  be  written  as 

6U  =  6w  F  (3.68) 

s 

where  6 w  is  the  virtual  displacement  and  F^  is  the  force  in  the 

spring.  Since  the  spring  is  linear  F  is  written  as 

s 

F  =  wkr/i  (3.69) 

s 

for  a  one  radian  sector  where  k  is  the  spring  stiffness  and  £  is 
the  length  and  r  is  the  radius  of  the  point  at  which  the  boundary 
element  is  attached.  If  the  spring  is  inclined  6°  to  the  vertical 
axis,  the  displacement  w  is  written  in  terms  of  the  horizontal  and 
vertical  components  as 


w 


u  sin0  +  v  cos0 


(3.70) 


' 
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Substituting  the  right  hand  sides  of  Eqs.  3.70  and  3.69 
into  Eq.  3.68  the  latter  is  written  in  matrix  form  as 


sin0 

u 

6u  =  <  6u  Sv  >  . 

■  <  sin0  cos0  >  Kr 

COS0 

V 

The  contribution  to  the  variation  of  the  total  internal 
virtual  work  of  the  structure  due  to  the  boundary  element  is  therefore 


u 


6U  =  <  6u  6v  >  [  ] 


(3.72) 


v 

where  [  K^e  ]  is  the  stiffness  matrix  of  the  boundary  element  and  is 
written  as 


[K 


sin20 
cosQ  sin0 


cos0  sin0 
cos20 


K  r 


(3.73) 


This  stiffness  matrix  is  added  to  the  appropriate  stiffness  coefficents 
of  the  point  (node)  at  which  the  boundary  element  is  attached. 

The  reactions  are  calculated  by  multiplying  [  ]  by  the 

actual  displacements  obtained  at  any  stage  in  the  solution.  The 
reactions  are  necessary  in  the  equilibrium  check.  They  are  added  to 
the  equilibrating  loads  obtained  in  section  3.4.3  as  a  method  of 
applying  the  condition  stated  by  Eq.  3.4b  to  virtual  displacement 
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(a)  Linear 


(b)  Quadratic 


(c)  Cubic 


Fig.  3.2  The  Serendipity  Finite  Element  Family 


Fig.  3.3  Normalized  Coordinates 
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(a)  Reinforcing  Layers 


(b)  The  Element  of  Length 

Fig.  3.4  The  Reinforcing  Element 


p,  |Pz 


(c)  Components 


Fig.  3.5  Surface  Tractions 


CHAPTER  FOUR 


PROGRAM  FE PARCS 5 


4 . 1  Introduction 

Program  FEPARCS5  is  a  finite  element  code  for  analysis  of 
axisymmetric  or  plane,  reinforced  and/or  prestressed  concrete  structure. 
In  this  program  it  is  assumed  that  displacements  are  small,  rotations 
are  negligible  and  strains  are  infinitesimal.  Although  the  program 
can  handle  linear  problems,  it  is  designed  for  problems  with  nonlinear 
material  response. 

In  this  chapter  a  general  description  of  the  program  is 
presented,  the  basic  solution  techniques  are  discussed  and  the  flow 
of  operations  for  the  program  is  briefly  described. 

4 . 2  General  Description 

Program  FEPARCS5  is  based  on  the  finite  element  model  formu¬ 
lated  in  Chapter  Three.  The  program  can  handle  combinations  of  linear, 
quadratic  or  cubic  isoparametric  elements.  The  solid  element 
formulated  in  Section  3.3.1  represents  the  concrete  continua. 

The  meridional  reinforcing  layers  and  prestressing  tendons  are 
represented  by  the  meridional  element  described  in  Section  3.3.2. 

Finally  the  element  formulated  in  Section  3.3.3  represents  the 
circumferential  reinforcing  layers  and  prestressing  tendons. 

The  constitutive  relation  proposed  in  Chapter  Two  for  three 
dimensional  (axisymmetric)  nonlinear  elastic  behavior  of  concrete  is 
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implemented  in  program  FEPARCS5  as  the  material  model  for  concrete. 

A  simple  one  dimensional  elastic  plastic  constitutive  relation  is 
adopted  for  the  reinforcing  layers  and  prestressing  tendons. 

The  program  can  handle  five  types  of  loads.  Each  type  is 
stored  in  a  separate  load  vector.  These  load  vectors  can  be  combined 
using  load  factors  specified  by  the  user  for  each  load  step  to  form  a 
load  increment  vector.  Dead  loads  can  be  a  combination  of  gravity 
loads,  hydrostatic  pressures,  and  concentrated  nodal  loads.  A 
separate  vector  for  live  concentrated  nodal  loads  is  provided.  A 
load  vector  is  provided  for  thermal  loads.  Normal  and  tangential 
surface  pressures  are  handled  by  two  separate  load  vectors.  Besides 
those  types  of  loads,  program  FEPARCS5  can  simulate  post  tensioning 
using  a  specially  prescribed  thermal  distribution.  This  procedure 
is  described  in  detail  in  Section  4.4.1. 

The  input  to  the  program  is  composed  of  control  parameters, 
material  properties,  nodal  geometry,  boundary  conditions,  solid 
element  information  and  topology,  reinforcing  and  prestressing 
layer  information,  concentrated  nodal  loads,  normal  and  tangential 
surface  pressure  nodal  intensity  distributions,  hydrostatic  pressure 
nodal  intensity  distribution  and  a  thermal  distribution.  The  speci¬ 
fication  of  each  load  combination  is  done  separately  for  each  load 
step . 

The  output  is  composed  of  nodal  displacements  and  incre¬ 
ments  of  displacements  in  the  global  coordinates,  local  coordinate 
stress  components  (in  the  and  ^ 1  directions  of  Fig.  3.3)  for 
solid  elements  and  stresses  and  strains  of  the  reinforcing  layers 


and  prestressing  tendons. 
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The  program  uses  numerical  integration  for  the  evaluation  of 
the  different  element  relations  as  well  as  the  loads  wherever  possible. 
A  number  of  Gaussian  integration  rules  have  been  provided  ranging 
from  one  point  rule  for  a  linear  displacement  four  node  element  to 
a  three  by  seven  two  dimensional  rule  for  higher  order  elements. 

The  program  can  use  a  tangential  stiffness  approach  to  the 
solution  or  alternatively  the  initial  load  method.  Both  these 
methods  will  be  described  in  Section  4.3.1.  The  equation  solving 
is  done  by  a  skyline  equation  solving  in-core  package  (see  Elwi, 

1977  and  Wilson  and  Bathe,  1975).  Dynamic  dimensioning  is  imple¬ 
mented  through  a  data  managing  package  described  by  Elwi  (1977) . 

Element  shape  functions  and  derivatives  evaluated  at  the  integration 
points  as  well  as  all  integration  point  information  such  as  stresses, 
strains,  and  material  properties  are  stored  on  sequential  files. 

In  this  manner  the  size  of  storage  required  by  the  program  is  kept 
to  a  minimum. 


4 . 3  Basic  Solution  Techniques 

4.3.1  The  Numerical  Solution  Strategy 


In  Chapter  Two  a  three  dimensional  constitutive  relation 
for  concrete  was  proposed.  This  relation  can  be  summarized  as 


dQ 


F  (d£ 


(4.1) 


in  which  F  represents  a  functional  relationship,  (Elwi  and  Murray, 
1979).  Eq.  4.1  implies  that  the  stress  increment  is  dependent  on 
the  previous  stress  history.  In  Chapter  Three  a  finite  element 


< 
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model  was  formulated  for  axisymmetric  reinforced  concrete  structures. 
This  model  is  a  displacement  model.  Therefore,  it  satisfies  kine¬ 
matic  compatibility  everywhere  and  it  approximately  satisfies 
equilibrium  only  on  a  global  level.  The  incremental  variational 
formulation  included  in  Section  3.2  leads  to  the  following  set  of 
equations . 

[  K  ]  {  Ar  }  =  {  AR  }  -  {  AQ  }  (4.2) 

which  is  identical  to  Eq.  3.16.  This  equation  is  considered  a  piece- 
wise  linearization  of  the  nonlinear  structural  response. 

Since  {Ar}  appearing  in  Eq.  4.2  is  of  finite  predetermin¬ 
ed  magnitude,  the  result  of  solving  the  set  of  equations  is  an 
increment  of  displacement.  The  increment  of  displacement  yields 
an  increment  of  strain  which  through  the  constitutive  matrix  yields 
an  increment  of  stress.  The  difference  between  the  applied  loads 
and  the  equilibrating  loads,  equivalent  to  the  stress  state  which 
satisfies  the  constitutive  law,  is  called  the  unbalanced  load.  One 
way  to  arrive  at  the  state  of  stress  which  is  kinematically  compatible 
and  at  the  same  time  satisfies  equilibrium,  is  to  eliminate  the 
unbalanced  load  through  an  iterative  scheme.  Iterative  schemes  of 
the  Newtonian  type  can  be  divided  into  two  main  categories.  These 
are  the  tangential  stiffness  method  (Argyris,  et  al.,  1974)  which 
is  sometimes  known  as  the  Newton  Raphson  Method  (Zienkiewicz ,  1971) 
and  the  initial  load  method  (Argyris,  et  al.,  1974)  which  is  sometimes 
known  as  the  modified  Newton  Raphson  method  (Zienkiewicz,  1971) . 

In  the  tangential  stiffness  method,  a  new  stiffness  matrix 
is  evaluated  at  the  beginning  of  each  load  increment,  based  on  the 
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current  material  properties.  This  method  is  powerful  and  comes  near 
to  simulating  the  actual  load  process.  The  main  disadvantage  of 
this  method  appears  when  strain  softening  behavior  is  exhibited.  In 
this  case,  the  tangential  stiffness  matrix  may  lose  its  positive 
definite  character  or  at  least  become  ill  conditioned.  In  this 
study,  however,  it  has  been  assumed  that  the  amount  of  steel  rein¬ 
forcing  and  prestressing  tendons  is  enough  to  retain  the  positive 
definiteness  of  the  problem.  In  addition,  the  stiffness  matrix 
may  become  nonsymmetric  when  nonhomogeneous  materials  are  present 
or  where  nonassociated  flow  rules  are  used  (Argyris,  et  al.,  1974). 
In  this  case,  the  computational  effort  associated  with  re-evaluating 
the  stiffness  matrix  may  become  prohibitive. 

In  the  initial  load  method  a  positive  definite  stiffness 
matrix  (usually  the  initial  stiffness  matrix)  is  retained  and  the 
softening  is  simulated  by  altering  the  unbalanced  load  vector 
(Argyris,  et  al . ,  1974).  Although  this  method  avoids  some  of  the 
problems  associated  with  the  tangential  stiffness  method,  a  much 
larger  number  of  iterations  may  be  required  to  satisfy  equilibrium 
particularly  when  approaching  the  ultimate  load. 

The  two  methods  described  above  can  be  illustrated  schema¬ 
tically  as  shown  in  Figures  4.1a  and  4.1b.  Many  improvements  can 
be  introduced  for  both  methods.  For  example,  over-relaxation  may  be 
used  to  improve  the  convergence  rate  of  the  initial  load  method. 
Under-relaxation  and  numerical  damping  can  be  used  to  enlarge  the 
convergence  domain  of  the  tangential  stiffness  method  (Almroth, 

Stren  and  Brogan,  1979  and  Fellipa,  1974  and  1976).  Re-evaluation 
of  the  stiffness  matrix  every  few  iterates  within  a  load  step  can 
be  carried  out  in  association  with  the  tangential  stiffness  method. 
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In  program  FEPARCS5  both  the  tangential  stiffness  approach 
with  the  capability  of  re-evaluating  the  stiffness  matrix  every  few 
iterates  and  the  initial  load  approach  have  been  implemented. 
Provisions ;  have  been  made  for  a  relaxation  factor.  Convergence  of 
the  iterative  scheme  is  based  on  both  the  displacements  and  the 
loads.  The  test  for  convergence  is  carried  out  as  follows 


<  X 

—  r 


(4.3a) 


<  X 

—  R 


(4.3b) 


where  X^_  and  X^  are  the  user  specified  tolerances  on  the  displace¬ 
ments  and  loads  respectively,  Ar1  is  the  increment  of  displacement 

th 

vector  obtained  in  the  l  iteration,  r  is  the  current  displacement 
vector,  Aq1  is  the  unbalanced  load  at  the  end  of  the  i*'*1  iteration 
and  R  is  the  total  load  vector.  The  symbol  | |  | |  denotes  the  Eucleadian 

norm. 


4.3.2  The  Subincrement  Method 

The  piecewise  linearization  process  described  in  Section 
4.3.1  requires  that  the  load  increment  be  reasonably  small.  When 
the  constitutive  law  is  path  dependent  the  resulting  strain  increments 
may  still  be  large  enough  to  cause  drift  of  the  solution, particularly 
when  strain  softening  behavior  is  exhibited.  The  ultimate  strength 
parameters  defined  in  Section  2.3.4  lie  on  the  intersection  of  a 
stress  path  with  the  ultimate  surfaces  described  in  Section  2.3.7. 
Complicated  ultimate  surfaces  may  yield  false  parameters  or  roots 
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outside  the  defined  range  when  a  stress  point,  which  does  not  yet 
satisfy  the  constitutive  relation,  falls  outside  the  surface. 

The  subincrement  method  avoids  these  problems  by  dividing 
the  strain  increments  into  a  number  of  smaller  sub increments .  In 
this  manner,  the  stress  point  is  changed  gradually  allowing  close 
simulation  of  the  behavior  and  the  constitutive  law  remains  stable 
even  when  the  first  derivative  of  the  stress  strain  curve  is  not 
linear.  Many  investigators  have  used  this  method  successfully  such 
as  Stricklin,  Haisler  and  von  Reissmann  (1971),  Bushnell  (1973) 
and  Murray,  et  al.  (1978) . 

While  the  subincrement  method  ensures  convergence  to  the 
right  answer,  it  is  expensive  because  the  material  handling 
routines  consume  a  substantial  portion  of  the  execution  time  in 
nonlinear  programs.  Therefore,  choosing  the  size  of  the  subincrements 
must  be  done  carefully.  Murray,  et  al.  (1978)  have  used  a  scheme 
based  on  the  size  of  the  strain  increment  as  compared  to  the 
strain  range.  In  FEPARCS5 ,  however,  the  number  of  subincrements 
is  set  by  the  user  for  each  load  step.  Thus,  the  economy  of  the 
solution  is  closely  controlled. 

4.3.3  Implementation  of  the  Constitutive  Relation 

The  implementation  of  the  constitutive  relation  proposed  in 
Chapter  Two  is  divided  into  two  main  parts.  In  the  first  part  the 
ultimate  surfaces  proposed  in  Section  2.3.7  are  prepared  by 
evaluating  the  constants  of  Eqs .  2.43a  and  2.43b  defined  in 
Appendix  A.  This  part  is  carried  out  in  the  problem  preparation 


phase  of  FEPARCS5. 


. 


72 


The  second  part  of  the  implementation  of  the  constitutive 
relation  is  related  to  active  usage  during  the  solution  phase.  In 
this  part,  the  increments  of  strain  are  evaluated,  the  equivalent 
uniaxial  strains  are  updated  and  the  corresponding  material  properties 
and  stresses  are  calculated  at  all  integration  points.  Rather  than 
using  the  increments  of  displacement  resulting  from  the  solution 
of  Eq.  4.2  in  each  iterate,  the  increments  are  added  to  form  a  total 
increment  of  displacement  vector  starting  at  the  beginning  of  each 
load  step.  This  vector  is  then  used  to  trace  the  variation  of  the 
stress  point  from  the  last  converged  position.  In  this  manner 
drift  of  the  solution  is  minimized.  In  the  following  the  mechanics 
of  this  phase  will  be  described  in  detailed  steps. 

1.  The  strain  increment  is  calculated  using  the  shape 
functions  and  derivatives  as 

{  Ae  }  =  [  B  ]  {  Aq  }  (4.4) 

2.  The  strain  subincrement  is  calculated  as 

<  Ae  >  =  <  Ae  >/N  (4.5) 

s 


where  N  is  the  number  of  subincrements. 

3.  The  strain  subincrement  is  then  transformed  into  the 
integration  point  local  coordinates  (chosen  in  Section  2.3.3  as 
the  fixed  orthotropy  axes)  as 


<Ae>* 


<  Ae  >  -*■  local  coordinate  system 

s 


(4.6) 
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The  program  then  loops  over  the  number  of  subincrements  to  perform 
the  next  steps . 

4.  The  local  stresses  are  updated  as 

(a  }J+1  =  {  a  +  [  c  ]  {  Ae}  £  (4.7) 


where  "i"  indicates  the  number  of  subincrement  which  varies  from  "o" 
to  "N-l"  with  "o"  denoting  the  state  of  stress  at  the  end  of  the 
previous  load  increment. 

5.  The  equivalent  uniaxial  strains  are  updated  as 


<  £  > 
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i+i 


=  <  £ 
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-  aJ)/E1> 


(4.8) 


6.  The  strength  and  deformation  parameters  are  calculated 


for  the  stress  point  <  0 ^  > 


i+i 


and  the  ultimate  surfaces 
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stress  ultimate  surface 


(4.9a) 
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■+  strain  ultimate  surface 


(4.9b) 


7.  The  stress  point  is  now  updated  using  the  stress  equi¬ 
valent  uniaxial  strain  curves  described  in  Section  2.3.4,  the 
material  properties  are  updated  using  the  functions  described  in 
Sections  2.3.5  and  2.3.6  and  the  constitutive  matrix  is  updated 
using  the  form  described  in  Section  2.3.2  as 

i+i 

<<J„  >  =  <fi  (£  ,  £  ,  a  ,  E  )  >  (4.10a) 

£  u  c  c  o 

<E>1+1  =  <fo  (£  ,  £  ,  G  ,  E  )>  (4.10b) 

^  u  c  c  o 

i+ 1 

< V  >  =  <f, (£  ,  e  ,  v  )> 

3  U  C  O 


(4.10c) 
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[C]  =  [f „(<E>1+1,  <V>1+1)]  (4 . lOd) 

in  which  f  denotes  functions. 

8.  At  the  end  of  the  loop  over  the  subincrements  the  program 
transforms  the  local  stresses  into  the  global  coordinate  system  and 
obtains  principal  stresses. 

The  eight  steps  described  above  are  repeated  at  each  inte¬ 
gration  point  for  all  elements.  It  must  be  noted  that  the  constitu¬ 
tive  matrix  has  been  formed  and  stored  in  the  local  coordinate 
system.  Therefore,  when  formulating  the  stiffness  matrix  of  the 
structure  care  must  be  taken  to  transform  the  constitutive  matrix 
into  the  global  coordinate  system. 

Initial  strains  caused  by  temperature,  if  present,  must  be 
removed  from  the  strain  increment  calculated  in  step  1  by  modifying 
Eq.  4.4  to  become 

{Ae}  =  [  B  ]  {Aq}  -  {  a}<0>{At}  (4.11) 

where  <  (J)  >  is  the  vector  of  shape  functions,  <  At  >  is  the  tempera¬ 
ture  increment  and  <  a  >  contains  the  thermal  expansion  coefficients 
which  are  assumed  not  to  be  affected  by  orthotropy. 

4 . 4  Special  Capabilities 
4.4.1  Initial  Stresses 

In  some  problems  where  a  self-equilibrating  stress  distribution 
exists,  it  becomes  necessary  to  evaluate  the  corresponding  strain 
field  before  the  solution  proceeds.  This  option  has  been  implemented 
in  FEPARCS5  in  the  initialization  phase.  Thus  instead  of  initializing 
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the  stresses  and  strains  to  zero  and  material  properties  to  the 
initial  state,  the  program  reads  the  stress  state  from  a  file  and 
calculates  a  corresponding  strain  field  that  satisfies  equilibrium 
but  may  not  satisfy  compatibility.  This  process  is  carried  out  in 
steps  similar  to  the  flow  chart  of  Fig.  2.16. 

4.4.2  Post  Tensioning 

Post  tensioning  is  a  common  prestressing  method  in  axisym- 
metric  structures.  The  method  consists  of  running  ducts  through 
the  structure  before  pouring  the  concrete.  When  the  concrete  is 
poured  and  has  attained  sufficient  strength,  the  prestressing  tendons 
are  inserted  in  the  ducts,  tensioned  and  anchored  to  the  concrete. 
During  this  process  the  structure  deforms  but  there  is  no  compatibility 
between  the  tendons  and  the  concrete.  The  next  step  is  to  grout 
the  ducts  to  protect  the  tendons  and  to  establish  compatibility 
of  deformation  between  tendons  and  concrete  under  all  subsequent 
loading  stages. 

Program  FEPARCS5  has  an  option  which  simulates  this  process 
very  successfully.  For  this  purpose  the  program  reads  a  special 
thermal  distribution  which  is  applied  only  to  the  prestressing  tendons 
causing  an  initial  strain  equal  in  magnitude  to  the  strain  caused  by 
the  prestressing  excluding  elastic  losses  (approximately) .  The 
program  then  assigns  each  tendon  a  level  of  prestressing  corresponding 
to  the  initial  strain.  Nodal  loads  are  then  calculated  from  the 
prestressing  elements  alone  using  Eqs .  3.57  and  3.58  of  Section  3.4.3. 

Since  at  that  stage  there  is  no  compatibility  between  tendons 
and  concrete,  the  program  evaluates  the  stiffness  matrix  of  the 
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structure  excluding  the  prestressing  elements.  The  loads  calculated 
from  the  prestressing  elements  are  then  subtracted  from  the  initialized 
load  vector  which  is  then  applied  to  the  structure.  At  the  end  of 
this  process  a  self-equilibrating  state  of  stress  exists  in  the 
structure  which  simulates  a  post  tensioning  procedure. 

This  process  is  considered  a  load  step,  the  end  results  of 
which  are  used  as  initial  conditions  for  the  next  load  step,  except 
that  the  displacements  and  loads  are  reassigned  to  the  state  before 
prestressing  once  convergence  is  obtained. 

4 . 5  Flow  of  Operations  for  Program  FE PARCS 5 

From  the  operational  point  of  view,  program  FEPARCS5  is 
divided  into  two  main  execution  stages,  namely  the  problem  prepara¬ 
tion  stage  and  the  solution  stage.  The  former  is  performed  only 
once  or  as  many  times  as  it  takes  to  check  the  data,  whereas  the 
latter  is  performed  for  each  load  step. 

As  outlined  in  Section  4.2  the  problem  preparation  stage 
reads  and  generates  the  structural  and  the  loading  data,  calculates 
and  stores  the  element  shape  functions  and  derivatives  at  all  Gaussian 
integration  points,  forms  the  sky  line  of  the  structure  stiffness 
matrix  and  initializes  the  stresses,  strain  and  material  properties 
at  all  integration  points.  Finally,  for  each  load  type  (dead  loads, 
live  concentrated  nodal  loads ,  thermal  loads ,  and  normal  and  tangen¬ 
tial  surface  pressure  loads)  a  separate  load  vector  is  formed.  If 
requested  an  initial  stiffness  matrix  is  formulated,  triangularized 


and  stored  to  be  used  in  the  initial  load  method. 
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Program  FEPARCS5  does  not  choose  the  size  of  the  load  step 
automatically.  This  is  done  by  the  user.  At  the  beginning  of  every 
load  step,  the  user  specifies  the  number  of  subincrements  and  whether 
the  tangential  stiffness  method  is  to  be  used  or  the  initial  load 
method.  If  the  former,  the  user  specifies  the  number  of  iterates 
after  which  the  stiffness  matrix  is  to  be  reevaluated.  Also  required 
are  the  tolerances  on  convergence,  the  relaxation  factor  and  finally 
load  factors  to  specify  the  load  combination  for  this  load  step. 

Having  read  the  user  load  step  specification  and  having 
formed  the  load  increment  accordingly,  the  program  proceeds  as 
follows . 

1.  The  program  formulates  and  triangularizes  the  stiffness 
matrix  in  the  case  of  the  tangential  stiffness  method  or  reads  the 
initial  triangularized  stiffness  matrix  from  out  of  core  in  the 
case  of  the  initial  load  method. 

2.  The  program  then,  solves  for  an  increment  of  displacement, 
updates  the  total  displacement  vector  and  updates  the  total  increment 
of  displacement  vector. 

3.  The  program  updates  the  stresses  and  material  properties 
according  to  the  scheme  described  in  Section  4.3.3. 

4.  If  the  problem  is  linear,  the  program  outputs  the 
results.  If  the  problem  is  nonlinear  the  program  calculates  the 
equilibrating  loads  using  the  equations  of  Section  3.4.3,  obtains 
the  unbalanced  load  vector  and  checks  for  convergence  as  described 
in  Section  4.3.1. 

5.  If  convergence  has  been  obtained  the  results  are  output 
and  the  load  step  is  considered  ended.  If  convergence  has  not  been 
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obtained,  and  the  number  of  iterates  has  not  exceeded  the  specified 
maximum,  the  program  uses  the  unbalanced  load  vector  to  obtain  a 
further  increment  of  displacement.  For  this  purpose  the  program 
branches  back  to  either  step  1,  if  re-evaluation  of  the  stiffness 
matrix  is  required  or  step  2  if  this  is  not  required. 

If  numerical  problems,  such  as  an  ill  conditioned  stiffness 
matrix,  oscillatory  convergence,  or  exceeding  the  maximum  number  of 
iterates,  occur,  the  program  automatically  stops  and  prints  out  the 
current  state  of  stress  and  displacement  for  the  user's  consideration. 

Fig.  4.2  and  4.3  show  the  flow  charts  of  the  operations  for 
the  two  stages  described  above. 


» 


Fig.  4.1  Iterative  Schemes 


Fig.  4.2  Flow  Chart  of  Problem  Preparation  Stage  of  FEPARCS5 
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CHAPTER  FIVE 


PRELIMINARY  INVESTIGATION 


5 . 1  Introduction 

A  constitutive  relation  for  three  dimensional  (axisymmetric) 
behavior  of  concrete  has  been  proposed  in  Chapter  Two.  In  Chapter 
Three  a  finite  element  model  for  incremental  analysis  of  axisymmetric 
reinforced  and/or  prestressed  concrete  structures  has  been  formulated. 
Chapter  Four  describes  program  FEPARCS5  in  which  the  features 
proposed  and  formulated  in  Chapters  Two  and  Three  have  been  imple¬ 
mented.  In  this  Chapter  a  preliminary  investigation  of  the  capabili¬ 
ties  of  program  FEPARCS5  and  its  ability  to  characterize  the  axisym¬ 
metric  behavior  of  concrete  is  conducted. 

A  series  of  tests  on  prestressed  concrete  wall  segments 
were  carried  out  at  the  I.F.  Morrison  Laboratory,  University  of 
Alberta  by  Dr.  J.G.  MacGregor,  Dr.  S.  Simmonds,  Dr.  D.W.  Murray  and 
Dr.  S.  Rizkallah  during  the  year  1978  (Chitnuyanondh ,  et  al.,  1979). 
These  experiments  formed  a  part  of  the  research  program  into  the 
effects  of  overpressure  on  secondary  containment  structures  mentioned 
in  Section  1.1  In  this  Chapter  a  simulation  of  some  of  these  tests 
is  carried  out  using  program  FEPARCS5 .  The  results  are  compared 
with  the  experimental  results  and  with  the  results  of  an  elastic 
plastic  analysis  (Chitnuyanondh,  et  al . ,  1979).  Only  two  segments 
out  of  fourteen  will  be  compared.  These  are  segments  no.  1  and  no.  3. 
Segments  no.  1  and  no.  3  are  identical  in  details  but  have 
been  tested  under  different  load  paths. 
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5 . 2  Wall  Segment  Test  Description 

The  wall  segment  specimens  were  31.5  x  31.5  inches  and  10.0 
inches  thick.  The  reinforcing  consisted  of  two  #3  @3  inches  meshes 
placed  one  half  inch  from  each  face.  The  prestressing  tendons 
consisted  of  four  7  wire  tendons  in  one  direction  and  three  6  wire 
tendons  in  the  other  direction.  The  tendons  were  placed  in  ducts 
1.07  inches  O.D.,  tensioned  and  anchored  at  the  ends.  Load  cells 
were  placed  between  the  anchor  heads  and  the  bearing  plates.  Figures 
5.1  to  5.3  show  the  details  of  the  wall  segments. 

The  primary  instrumentation  consisted  of  electric  strain 
gages,  placed  on  the  concrete  on  each  face  parallel  to  the  main  direc¬ 
tions  of  prestressing.  A  number  of  Dsmec  points  were  also  installed 
on  both  faces. 

Loading  was  applied  through  two  independent  perpendicular 
systems  attached  to  both  the  steel  reinforcing  and  the  prestressing 
tendons.  Segment  1  was  tested  under  biaxial  tension  with  a  2:1 
loading  ratio.  Segment  3  was  also  tested  under  biaxial  tension  but 
the  loading  ratio  was  1:1  up  to  375.0  kips  at  which  point  the  load  in 
direction  2  was  kept  constant  while  the  load  in  direction  1  was 
increased.  Direction  1  was  the  direction  parallel  to  the  7-wire 
tendons,  while  direction  2  was  that  parallel  to  the  6-wire  tendons. 

The  prestressing  level  in  the  tendons  in  direction  1  was  134.6  psi 
and  in  direction  2  it  was  124.2  ksi  (losses  included).  The  prestress¬ 
ing  system  was  identical  in  both  segments  1  and  3. 
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5 . 3  FEPARCS5  Modelling  of  the  Wall  Segments 

Wall  segments  no.  1  and  no.  3  have  been  modelled  for  analysis 
by  program  FEPARCS5  as  parts  of  a  cylinder  31.5  inches  high  and  10.5 
inches  thick  having  an  110.5  inches  internal  radius  as  shown  in 
Fig.  5.4.  The  finite  element  model  consists  of  one  12  node  cubic 
element  resting  on  roller  boundary  conditions  as  shown  in  Fig.  5.5. 

The  steel  rebars  and  prestressing  tendons  are  smeared  in  meridional 
and  circumferential  layers  as  shown  in  Fig.  5.6. 

The  material  properties  are  based  on  data  published  by 
Chitnuyanondh ,  et  al .  ,  (1979).  Table  5.1  contains  the  material 
properties  required  for  the  proposed  constitutive  relation.  Table  5.2 
describes  the  stress  strain  curves  for  the  prestressing  tendons  and 
the  steel  rebars. 

Surface  pressure  is  applied  on  the  internal  face  of  the  cylinder 
in  order  to  induce  tensile  stresses  in  the  circumferential  direction 
which,  henceforth,  is  called  direction  2.  A  consistent  set  of 
concentrated  nodal  loads  is  applied  to  the  top  edge  in  order  to  induce 
tensile  stresses  in  the  meridional  direction  which,  henceforth,  is 
called  direction  1.  The  prestressing  is  applied  using  the  post¬ 
tensioning  option  described  in  Section  4.4.2.  The  tangent  stiffness 
approach  is  used  throughout  the  analysis  of  both  segments  as  a 
solution  strategy. 

5.4  Comparison  of  FEPARCS5  Analysis  with  the  BOSOR5  Theoretical 
Analysis  and  with  the  Experimental  Analysis 

As  mentioned  in  Section  5.2,  segment  no.  1  has  been  loaded 
with  a  ratio  of  2:1.  Fig.  5.7  shows  the  load-strain  response  of 


•• 

* 

■ 


85 


this  segment.  On  this  figure,  the  response  of  FEPARCS5  in  direction  1 
agrees  well  with  that  of  B0S0R5  and  with  the  experimental  results. 

In  Direction  2  FEPARCS5  and  B0S0R5  agree  well  but  both  predict 
higher  cracking  load  and  subsequently  a  response  stiffer  than  that 
of  the  experiment. 

Segment  no.  3  has  been  loaded  with  a  1:1  ratio  up  to  375.0 
kips,  at  which  point  the  loading  in  direction  2  has  been  maintained 
constant,  while  the  loading  in  direction  1  has  been  continued. 

Fig.  5.8  shows  the  load-strain  response  of  this  segment.  In  this 
case  FEPARCS5  shows  a  more  flexible  response  in  direction  1  compared 
to  both  B0S0R5  and  the  experiment.  Direction  2  shows  good  agreement 
between  all  three  responses.  Table  5.3  contains  the  loads  at  which 
cracking,  rebar  yield  and  tendon  yield  have  occurred. 

5 . 5  The  Stress  Path 

The  stress  path  observed  for  segment  no.  3  gives  strong 
evidence  supporting  the  assumption  of  Section  2.3.8  that  in  biaxial 
and  triaxial  tension  cracking  in  one  direction  does  not  precipitate 
cracking  in  the  perpendicular  direction.  This  stress  path  is  shown 
in  Fig.  5.9.  Fig.  5.9  indicates  that  although  direction  2  reaches 
the  maximum  tensile  strength  and  starts  to  descend,  direction  1 
continues  to  accept  increasing  stresses  until  it  assumes  the  maximum 
uniaxial  tensile  strength.  At  this  point  cracking  normal  to  direction 
2  occurs  and  the  stress  in  this  direction  starts  to  decrease.  The 
stress  paths  for  segment  no.  3  obtained  from  the  experiment  and 
B0S0R5  analysis  (Chitnuyanondh)  are  shown  on  Fig.  5.9  together  with 
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that  of  FEPARCS5.  Agreement  between  all  three  paths  is  fair  as 
far  as  values  are  concerned.  The  pattern  discussed  above  however  is 
very  distinct  in  all  three  paths. 
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Ultimate  Strength 
Surface  Parameters 

Equivalent  Uniaxial 

Strain  Surface  Parameters 

Elastic  Moduli 

f  * 

-5090.0 

e 

0.002 

E,  ,  * 

4.1  x  106 ’ 

cu 

cu 

d  l 

a 

1.2 

a 

1.3 

f_  _  * 

4.1  x  106 

c 

c 

0  2 

at 

0.034 

at 

0.095 

^0  3 

4.1  x  106 

13.75 

Si 

22.5 

G0  12  * 

1.7  x  106 

Pi 

0.0 

Pi 

0.0 

V0  1 

0.2 

s2 

3.75 

S2 
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^0  2 

0.2 

P2 

0.0 

P2 

0.0 

V0  3 

0.2 

(a)  Wall  Segment  No.  1 


Ultimate  Strength 
Surface  Parameters 

Equivalent  Uniaxial 

Strain  Surface  Parameters 

Elastic  Moduli 

f  * 

-5690.0 

£ 

0.002 

En  i 

4.3  x  106 

cu 

cu 

U  I 

a 

c 

1.2 

a 

1.3 

Eo  2 

2.2  x  106 

at 

0.034 

a 

0.12 

E0  3 

2.2  x  106 

Si 

13.75 

Si 

22.5 

^  12 

1.2  x  106 

Pi 

0.0 

Pi 

0.0 

V0  1 

0.2 

C  2 

3.75 

C2 

22.5 

^0  2 

0.2 

P2 

0.0 

p2 

0.0 

V0  3 

0.2 

(b)  Wall  Segment  No.  3 


Table  5.1  Wall  Segment  Material  Parameters  (Concrete) 


#3  Rebars 

Prestressing  Tendons 

0 

£ 

0 

£ 

-3 

-  3 

(ksi) 

(10  ) 

(ksi) 

(10  3) 

0 

0 

0 

0 

58.20 

2.04 

205.0 

6.97 

69.60 

40.00 

228.0 

8.40 

— 

— 

237.0 

10.0 

240.0 

12.0 

250.0 

20.0 

251.2 

41.0 

Table  5.2  Wall  Segment  Material  Parameters 
(Steel) 
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Wall  Segment  No.  1 

Wall  Segment  No. 

3 

FEPARCS5 

B0S0R5 

TEST 

FEPARCS5 

B0S0R5 

TEST 

Cracking  of 
Direction  1 

310 

295 

320 

300 

263 

330 

Cracking  of 
Direction  2 

430 

350 

359 

225 

200 

207 

Vertical  Rebar 

Yield 

490 

475 

465 

500 

530 

530 

Horizontal  Rebar 
Yield 

— 

— 

— 

350 

350 

355 

Vertical  Tendon 
Yield 

510 

495 

500 

500 

— 

— 

Horizontal 

Tendon  Yield 

— 

— 

— 

375 

425 

375 

Table  5.3  Limit  State  Loads  (all  loads  shown  are 
those  in  direction  1  in  kips) 


90 


Fig.  5.1  Elevation  of  Wall  Segment  Specimen 
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Fig.  5.2  Typical  Section  A-A  Through  Wall  Segment  Specimen 


Fig.  5.3  Typical  Section  B-B  Through  Wall  Segment  Specimen 
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Direction  1 


Fig.  5.4  Axisymmetric  Modelling  of  Wall 
Segment  Specimen 


Fig.  5.5  Finite  Element 

Modelling  of  Wall 
Segment  Specimen 


Fig.  5.6  Reinforcing  and 

Prestressing  Layers  of 
Wall  Segment  Finite 
Element  Model 
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Fig.  5.7  Load  Strain  Comparison  for  Segment  No. 
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Fig.  5.8  Load  Strain  Comparison  for  Segment  No.  3 
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Fig.  5.9  Stress  Path  Comparison  for  Segment  No.  3 
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CHAPTER  SIX 


ANALYSIS  OF  A  CONTAINMENT  STRUCTURE 


6.1  Introduction 


In  Chapter  Five  a  preliminary  investigation  of  the  capabili¬ 
ties  of  program  FEPARCS5  has  been  undertaken.  The  next  step  is  to 
use  the  program  to  analyse  a  complicated  structure  which  makes  more 
rigorous  demands  on  the  different  features  of  the  program.  In  this 
chapter  a  finite  element  model  of  the  test  structure  mentioned  in 
Section  1.1  is  described,  an  analysis  of  the  model  using  program 
FEPARCS5  is  presented,  and  a  comparison  of  the  results  both  with 
those  obtained  from  the  actual  test  and  with  those  of  an  elastic 
plastic  theoretical  analysis  of  a  similar  model  carried  out  by 
Murray,  et  al.  (1978)  using  a  modified  form  of  program  B0S0R5  is 
discussed. 

6 . 2  Description  of  Test  Structure  and  Procedure 

The  test  structure,  shown  in  Figures  6.1  to  6.4,  was  composed 
of  a  3' -6"  thick  base,  a  5"  thick  -  10'  high  cylindrical  wall  which 
had  an  internal  radius  of  4 ’-10",  a  ring  beam  and  a  4"  thick 
spherical  dome  which  had  an  internal  radius  of  9' -8".  A  smooth 
transition  between  the  dome  and  ring  beam  was  obtained  by  gradually 
thickening  the  dome  at  the  springing  line.  Four  buttresses  were 
placed  at  90°  intervals  around  the  cylinder  to  provide  anchorage 
for  the  circumferential  post-tensioning  strands.  The  dome  prestressing 
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tendons  were  placed  on  two  orthogonal  geodetic  nets  as  shown  in  the 
upper  right  hand  quadrant  of  Fig.  6.2.  The  dome  reinforcing  was 
placed  on  meridional  and  circumferential  lines  as  shown  in  the  lower 
right  hand  quadrant  of  Fig.  6.2. 

The  levels  of  prestressing  and  the  thicknesses  have  been 
designed  to  yield  a  cracking  sequence  similar  to  that  expected  in 
the  Gentilly-2  containment  structure,  (Murray  et  al.  1978). 

However,  the  wall  has  been  designed  to  have  lower  cracking  strength 
relative  to  that  in  the  dome  than  would  be  expected  in  the  Gentilly-2 
structure.  In  this  manner  more  information  can  be  obtained  about 
the  behavior  of  all  components  of  the  test  structure  before  final 
failure  occurred. 

Testing  of  the  structure  was  carried  out  by  water  pressure 
applied  through  .a  thin  elastic  membrane  which  prevented  leakage 
through  the  wall  of  the  structure.  Five  preliminary  tests,  in  which 
the  pressures  were  kept  within  20.0  psi,  were  carried  out  to  check 
instrumentation.  Subsequently,  two  actual  tests  were  conducted.  In 
the  first  test,  pressure  was  raised  up  to  80.0  psi  which  was  well 
above  cracking  pressure.  The  structure  was  then  unloaded  since 
substantial  leaking  of  water  through  a  damaged  membrane  made  it 
impossible  to  maintain  pressure.  A  new  liner  was  installed  and  the 
final  test  was  carried  out  up  to  failure  at  159.9  psi. 

First  visible  cracks  were  observed  at  40.0  psi  in  the  cylinder 
wall  along  meridional  and  circumferential  lines.  At  140.0  psi 
splitting  cracks  at  the  anchorage  of  a  circumferential  tendon  ap¬ 
peared  on  the  southwest  buttress.  At  158.0  psi  failure  occurred  at 
this  location  with  concrete  spalling  off.  At  159.0  psi  three 
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horizontal  tendons  in  the  wall  either  fractured  or  lost  their 
anchorage  and  the  adjacent  rebars  fractured.  The  liner  then 
ruptured  releasing  all  pressure  and  ending  the  test. 

Electronic  readings  were  taken  at  intervals  of  internal 
pressure  of  5.0  psi.  There  were  207  electric  strain  gages  placed  on 
the  steel  reinforcing  and  38  electric  strain  gages  placed  on  concrete 
face.  In  addition  74  Demec  gages  were  located  along  and  across  a 
meridional  line  midway  between  the  northeast  and  northwest  buttresses. 

The  information  on  the  test  contained  in  this  section  and 
the  next  sections  have  been  obtained  from  progress  reports  to  the 
Atomic  Energy  Control  Board  of  Canada,  since  final  reports  have  not 
been  published  up  to  the  time  of  writing. 

6 . 3  Finite  Element  Model  of  Test  Structure 
6.3.1  General  Description 

The  finite  element  mesh  for  the  structure  is  controlled  by 
the  locations  of  the  prestressing  strands,  the  shape  of  the  ring 
beam  and  the  anticipated  stress  gradients  at  the  junction  of  the 
ring  beam  and  the  cylinder  wall.'  Quadratic  displacement  element  are 
used  exclusively  to  construct  the  mesh.  These  elements  have  been 
found  to  give  greater  flexibility  in  modelling  the  geometry  of  the 
structure  than  do  linear  elements.  At  the  same  time  the  problem 
size  is  kept  reasonable.  Two  elements  are  used  through  the  thickness 
of  the  dome  and  the  cylinder  wall  to  enable  accurate  positioning  of 
the  prestressing  tendons  in  the  dome  and  to  resolve  the  problem  of 
intersection  of  vertical  and  dome  prestressing  tendons  inside  the 
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ring  beam.  Figures  6.5  to  6.8  show  the  mesh  and  the  prestressing 
and  reinforcing  layers  superimposed  on  the  solid  elements. 

It  has  been  found  by  Murray  et  al.  (1978) ,  that  the  base 
provides  complete  fixity  at  the  bottom  of  the  cylinder  wall. 

Therefore,  the  base  has  been  eliminated  from  the  finite  element  model 
and  in  its  place  a  set  of  boundary  elements  have  been  used  as  shown 
in  Figures  6.6b  and  6.6c.  Since  only  one  half  of  the  vertical  cross 
section  is  modelled,  the  end  of  the  dome  on  the  axis  of  symmetry  has 
been  provided  with  a  set  of  horizontal  boundary  elements  as  shown 
in  Fig.  6.6a  in  order  to  prevent  the  vertical  line  from  rotating  or 
translating  in  the  horizontal  direction. 

The  connectivity  of  the  elements  has  been  defined  so  that 
the  local  y-£  coordinates  be  as  shown  in  Fig.  6.9.  The  order  of 
the  Gaussian  integration  rule  is  2  x  2  for  all  solid  elements.  A 
two  point  rule  is  used  for  each  reinforcing  or  prestressing  layer 
within  a  solid  element.  Altogether,  the  problem  has  226  nodes,  14 
boundary  elements  and  55  solid  elements  on  which  a  total  of  63  longi¬ 
tudinal  reinforcing  layers,  57  hoop  reinforcing  layers,  29  longi¬ 
tudinal  prestressing  layers  and  28  circumferential  prestressing  layers 
are  imposed. 

6.3.2  Modelling  of  Dome  Prestressing 

In  program  FE PARCS 5  prestressing  layers  have  to  be  described 
as  meridional  or  circumferential  layers.  Therefore,  it  has  been 
necessary  to  transform  the  dome  prestressing  layers  from  the  actual 
orthogonal  geodetic  nets  into  equivalent  meridional  and  circumferen¬ 
tial  layers.  This  has  been  accomplished  by  deriving  weighting 
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functions  for  the  contribution  of  each  tendon  to  the  meridional 
and  circumferential  directions  at  the  points  of  intersection  of  the 
tendons.  The  weighting  functions  have  been  plotted  against  the  radii 
of  the  points  of  intersection  for  meridional  contributions  and 
circumferential  contributions.  These  plots  have  been  used  to  obtain 
the  average  contribution  of  meridional  tendons  per  unit  width  and 
circumferential  tendons  per  unit  length.  This  process  is  described 
in  detail  in  Appendix  C. 

6.3.3  Materials 

The  stress  strain  curves  for  the  different  kinds  of  steel 
rebars  are  shown  in  Figs.  6.10a  and  6.10b  (Murray,  et  al ,  1978). 

The  maximum  strain  allowed  in  rebars  is  4%.  The  stress  strain 
curves  for  the  prestressing  tendons  are  shown  in  Figs.  6.11a  and 
6.11b  (Murray,  et  al . ,  1978)  .  The  maximum  strain  allowed  in  the 
dome  tendons  (0.62"cj))  is  8%,  while  the  wall  prestressing  tendons 
(0.5"(J))  are  allowed  only  5%  strain. 

Two  different  types  of  concrete  are  used  in  the  body  of  the 
model.  Normal  cast  in  place  concrete  is  used  for  the  cylinder  wall 
up  to  an  elevation  of  5.0  ft.,  and  in  the  ring  beam  and  the  dome. 
Gunite  concrete  is  used  in  the  cylinder  wall  from  elevation  of  5.0  ft. 
to  10.0  ft.  The  properties  of  those  materials  are  shown  in  Table  6.1. 

No  thermal  analysis  of  this  model  is  required,  hence  concrete 
and  all  rebars  are  assigned  zero  thermal  properties.  However, 
since  the  prestressing  strains  in  the  tendons  are  assigned  through 
a  thermal  field,  the  tendons  are  assigned  thermal  properties  corres¬ 
ponding  to  the  required  levels  of  prestressing.  There  are  essentially 
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three  levels  of  prestressing:  137.7  ksi  for  the  circumferential 
prestressing  in  the  wall  and  the  ring  beam,  90.0  ksi  for  the  vertical 
prestressing  in  the  wall  and  113.1  ksi  for  the  dome  prestressing 
(losses  included).  The  first  two  levels  are  induced  in  0.5"(J)  tendons. 
Therefore,  two  different  material  types  are  defined  to  simulate  the 
0.5"(j)  tendons  with  two  different  thermal  expansion  coefficients.  The 
dome  prestressing  is  composed  of  0.62"(f)  tendons  and  these  are  defined 
to  have  a  different  thermal  expansion  coefficient. 

6.3.4  Loads 

This  model  is  subjected  to  a  variety  of  loads;  gravity  loads, 
tangential  friction  losses  in  the  dome,  hydrostatic  pressure, 
prestressing  loads  and  the  internal  pressure.  Gravity  loads  are 
simulated  as  described  in  Section  3.4.1.  Friction  losses  in  the  dome 
are  simulated  by  meridional  tangential  tractions  on  the  middle  line 
of  the  dome.  Hydrostatic  pressures  are  simulated  by  a  normal  pressure 
distribution  on  the  inside  face  of  the  model.  These  three  load 
vectors  are  added  to  one  dead  load  vector.  Prestressing  simulation 
has  been  carried  out  using  the  procedure  outlined  in  Section  4.4.2. 

A  normal  internal  pressure  of  1.0  psi  is  used  to  form  an  equivalent 
internal  pressure  load  vector  which  is  incremented  in  the  subsequent 
pressurization  of  the  structure. 

6 . 4  Analysis  of  Test  Structure 

The  analysis  of  the  finite  element  model  described  in  Section 
6.3  is  performed  in  several  stages.  In  the  problem  preparation  stage 
data  is  read  and  generated,  element  shape  functions  and  derivatives 
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are  evaluated,  column  heights  and  the  addressing  array  of  the  stiff¬ 
ness  matrix  are  formed,  stresses,  strains  and  material  properties 
are  initialized  at  all  integration  points,  and  all  basic  load 
vectors  are  formed.  In  the  second  stage  the  tendons  are  tensioned, 
equivalent  prestressing  loads  are  evaluated  and  applied  to  the 
structure  excluding  the  tendons.  In  the  third  stage  dead  loads 
composed  of  gravity,  friction  and  hydrostatic  loads  are  applied  to 
the  structure  as  one  load  increment.  In  the  subsequent  stages,  incre¬ 
ments  of  internal  pressure  are  added  to  the  structure  simulating 
the  pressurization  of  the  test  structure  which  ultimately  leads  to 
failure.  Failure  in  this  analysis  is  characterized  by  the  inability 
of  the  program  to  obtain  convergence  on  both  loads  and  displacements. 

Table  6.2  contains  a  breakdown  of  all  runs  on  the  Amdahl 
470/V7  computer  unit  at  the  University  of  Alberta,  Edmonton.  The 
number  of  iterates  required  for  each  load  step  and  the  CPU  time 
consumed  are  shown.  Table  6.2  also  shows  the  tolerances  on  the 
displacements  and  the  loads,  A ^  and  A^  respectively.  The  number  of 
subincrements  (NS)  vary  between  5  and  10.  The  maximum  number  of 
iterates  allowed  per  load  step  is  31  up  to  125.0  psi,  after  which  the 
limit  is  increased  to  70  iterates.  In  the  initial  stages  and  up 
to  40.0  psi  the  stiffness  matrix  is  evaluated  at  the  beginning  of 
each  load  step  and  re-evaluated  after  each  third  iterate.  In 
subsequent  stages,  the  material  properties  available  at  20.0  psi  are 
used  to  formulate,  triangularize  and  store  a  constant  stiffness 
matrix  which  is  used  for  all  load  increments  above  40.0  psi.  It  has 
been  found  that  this  stiffness  matrix  gives  better  convergence  than 


the  initial  stiffness  matrix. 
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6.5  Comparison  of  Results  with  the  Experimental  Results  and  the 
B0S0R5  Theoretical  Analysis 

Program  FEPARCS5  outputs  nodal  displacements,  stresses  at 
the  Gaussian  points  for  solid  elements,  and  stresses  and  strains  at 
the  Gaussian  points  for  the  rebar  and  prestressing  layers.  At  the 
time  of  writing,  the  experimental  results  and  the  B0S0R5  results 
available  consist  of  displacement  profiles  and  surface  strains 
obtained  from  Demec  readings.  Comparisons,  therefore,  will  be  re¬ 
stricted  to  displacement  history  and  profiles,  and  surface  strains. 

In  addition,  cracking  patterns  are  discussed. 

6.5.1  Pi splacements 

Deflection  patterns  determined  from  the  FEPARCS5  analysis 
are  reasonably  similar  to  those  of  B0S0R5  and  the  test.  Fig.  6.12 
shows  a  deflection  profile  at  120.0  psi.  It  can  be  seen  that  the 
results  of  FEPARCS5  on  the  cylinder  wall  represent  a  stiffer  behaviour 
than  for  the  other  two  sets  of  results.  On  the  dome,  however, 

FEPARCS5  and  B0S0R5  match  reasonably  up  to  a  point  two  feet  away  from 
the  crown.  At  this  point  FEPARCS5  maintains  the  displacement  gradient 
an  extra  short  distance,  thereby  approaching  the  test  results, 
while  B0S0R5  results  level  off  more  rapidly.  On  this  figure  it  is 
observed  that  the  displacement  gradient  of  the  test  structure  at 
120.0  psi  is  higher  than  either  that  of  BOSOR5  or  FEPARCS5. 

Fig.  6.13  shows  the  deflection-pressure  curves  at  two 
points  near  the  top  of  the  dome  and  the  middle  of  the  cylinder. 

These  are  the  points  of  maximum  deflection.  On  this  figure,  B0S0R5 
and  FEPARCS5  agree  well  up  to  80.0  psi  at  the  top  of  the  dome  and 
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up  to  110.0  psi  at  the  middle  of  the  cylinder.  After  80.0  psi 
FEPARCS5  deflection  response  at  the  top  of  the  dome  approaches  the 
test  results  while  B0S0R5  predicts  a  somewhat  stiffer  response.  This 
agrees  well  with  the  deflection  profiles  shown  in  Fig.  6.12. 

FEPARCS5  deflection  response  at  the  middle  of  the  wall  after  110.0 
psi  predicts  more  carrying  capacity  than  that  shown  by  B0S0R5.  The 
reason  must  be  that  B0S0R5  uses  a  better  descending  branch  representa¬ 
tion  of  the  tensile  stress  strain  curve  of  concrete.  The  difference, 
however,  is  within  2.5%. 

6.5.2  Surface  Strains 

As  mentioned  before  the  strains  plotted  from  the  test 
results  are  measured  from  D.emec  readings,  i.e.  are  measured  from 
displacements  over  a  finite  length.  Hence,  at  points,  where  there 
is  a  relatively  high  displacement  gradient,  a  discrepancy  must  show 
between  those  results  and  the  strain  predicted  by  FEPARCS5 .  The 
latter,  as  mentioned  in  Chapter  Three,  are  calculated  assuming 
infinitesmal  strain  and  negligible  rotations.  This  is  parti¬ 
cularly  noticeable  in  the  cylinder  wall  near  the  base  where  there 
is  a  point  of  inflection  and  rotations  are  not  negligible. 

Figs.  6.14  to  6.19  show  a  series  of  strain  histories  against 
internal  pressure  at  a  variety  of  points  on  the  outside  surface  of 
the  structure.  Fig.  6.14  shows  the  meridional  strains  at  a  point 
28"  above  the  base  of  the  cylinder.  The  response  of  FEPARCS5  agrees 
reasonably  with  that  of  B0S0R5  and  the  test  results  up  to  105  psi. 

At  this  point  the  response  of  FEPARCS5  shows  considerable  stiffening 
as  shown  by  the  curve  denoted  by  This  discrepancy  is  to  be 
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expected  since  this  curve  is  derived  from  an  infinitesmal  strain  field 
and  the  point  under  consideration  lies  very  near  the  point  of  in¬ 
flection  as  can  be  seen  in  Fig.  6.12.  The  curve  denoted  £^  shows 
finite  (Green's)  strains  computed  from  the  displacement  field  of 
FEPARCS5.  This  curve  reflects  the  influence  of  rotations  on  the 
strain  at  this  point. 

Fig.  6.15  shows  the  circumferential  strains  at  the  same  point. 
Agreement  here  is  reasonable  between  all  three  sets  of  results, 
although  the  stiff  behavior  predicted  by  FEPARCS5  in  the  cylinder 
wall  can  be  observed. 

Figs.  6.16  and  6.17  show  respectively  the  meridional  and 
circumferential  strains  at  a  point  65"  above  the  base  on  the  outer 
surface  of  the  cylinder.  In  the  neighborhood  of  this  point,  the 
displacement  gradient  is  very  small  (see  Fig.  6.12)  and  FEPARCS5 
predicts  a  meridional  strain  response  in  good  agreement  with  B0S0R5 . 
The  test  results,  however,  are  much  stiffer.  Confirmation  of  the 
test  results  may  be  obtained  once  the  electric  strain  gage  results 
are  published.  In  the  circumferential  direction  agreement  is  good 
between  all  three  sets  of  results,  although  the  slightly  stiffer 
wall  response  of  FEPARCS5  is  again  observed. 

Fig.  6.18  shows  the  meridional  strains  at  a  point  approxi¬ 
mately  8.5"  away  from  the  crown  on  the  outer  surface  of  the  dome. 

FE PARCS 5  response  indicated  by  exhibits  the  stiff  behavior 
characterizing  points  of  high  displacement  gradient.  The  curve 
indicated  £^  is  obtained  including  the  rotation  from  the  displacement 
field  of  FEPARCS5.  The  improvement  in  response  is  noticeable. 

Fig.  6.19  shows  the  circumferential  strain  at  the  same  point. 


106 


FEPARCS5  response  agrees  well  with  that  of  B0S0R5  and  the  test 
results . 

6.5.3  Cracking  Sequence 

Program  FEPARCS5  outputs  cracking  information  at  the  Gaussian 
points.  The  term  "cracking"  means  that  the  point  has  reached  its 
ultimate  strength  and  is  forced  to  deform  beyond  the  corresponding 
strain.  The  term  'vertical  crack'  is  used  to  describe  cracks  occurring 
due  to  circumferential  strains,  while  the  term  'horizontal  crack' 
is  used  to  describe  those  cracks  occurring  due  to  meridional  strains. 

Fig.  6.20  show  the  progress  of  vertical  cracks  as  predicted 
by  FEPARCS5 .  Vertical  cracks  occur  as  early  as  40.0  psi  at  the  crown 
of  the  dome  and  spread  slowly  outward  on  the  outer  surface  of  the 
dome.  Extensive  vertical  cracks  appear  on  the  entire  length  of  the 
cylinder  with  the  exception  of  the  upper  and  lower  two  feet  at  60.0 
psi.  As  pressure  increases  vertical  cracks  spread  up  and  down  on 
the  cylinder  but  never  reach  the  ring  beam  nor  the  base.  Throughout 
the  analysis  the  ring  beam  and  the  springing  line  of  the  dome  remain 
free  of  vertical  cracks.  B0S0R5  predicts  extensive  vertical  through 
cracks  in  the  cylinder  wall  at  62.0  psi  (Murray,  et  al . ,  1978). 

Vertical  cracks  in  the  dome  start  appearing  at  67.0  psi  and  spread 
slowly  in  the  B0S0R5  analysis.  However,  at  120.0  psi  the  pattern  of 
vertical  dome  cracks  predicted  by  B0S0R5  is  similar  to  that  of 
FEPARCS5.  Again  the  ring  beam  and  the  springing  region  of  the  dome 
remain  free  of  vertical  cracks  throughout  the  analysis  (Murray, 
et  al.,  1978).  In  the  test  the  first  visible  signs  of  vertical 
cracks  were  observed  at  40.0  psi  in  the  wall.  At  80.0  psi  these 
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cracks  became  as  extensive  as  has  been  predicted  by  both  B0S0R5  and 
FEPARCS5 .  The  ring  beam  and  the  springing  of  the  dome  remained  free 
of  vertical  cracks  to  the  end  of  the  test. 

Fig.  6.21  shows  the  progress  of  horizontal  cracking.  The 
onset  of  horizontal  cracks  is  predicted  by  FEPARCS5  at  50.0  psi  on 
the  inside  surface  of  the  ring  beam  and  the  springing  of  the  dome.  At 
60.0  psi  FEPARCS5  predicts  horizontal  cracks  at  upper  surface  of  the 
dome  at  and  near  the  crown  and  on  the  outside  surface  of  the  wall  two 
feet  below  the  ring  beam.  The  cracked  region  on  the  inside  surface 
of  the  ring  beam  continued  to  spread  up  and  down  slowly  penetrating 
the  ring  beam.  At  70.0  psi  extensive  horizontal  cracks  are  predicted 
in  the  dome  and  the  wall.  However,  the  upper  surface  of  the  ring 
beam  and  the  springing  of  the  dome  remain  crack  free.  At  120.0  psi 
horizontal  cracks  appear  everywhere  in  the  structure  except  in  the 
upper  surface  of  the  springing  of  the  dome  the  outer  portion  of  the 
ring  beam  and  the  outside  surface  of  the  wall  near  the  base.  B0S0R5 
predicts  the  first  horizontal  cracks  at  30.0  psi  on  the  inside 
surface  of  the  ring  beam.  The  outside  surface  of  the  dome  near  the 
crown  and  the  outside  surface  of  the  wall  under  the  ring  beam  are 
predicted  to  crack  horizontally  at  62  psi.  After  that  B0S0R5  shows 
a  pattern  of  horizontal  crack  progress  which  is  very  similar  to  that 
predicted  by  FEPARCS5 .  Test  results  indicate  similar  behavior  on 
the  outside  surface.  On  the  inside  surface,  no  information  is  avail¬ 


able  since  it  was  inaccessible. 
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6 . 6  Limit  States  for  Containment  Structures 

The  objective  of  the  application  of  the  technology  developed 
in  this  work  is  to  simulate  the  response  of  a  containment  structure 
to  internal  over  pressure.  Internal  pressure  in  nuclear  containment 
structures  may  occur  due  to  failure  of  the  primary  or  secondary 
cooling  systems  within  the  structure.  Under  such  hypothetical 
situations  a  series  of  limit  states  can  be  defined  to  indicate  the 
degree  of  damage  (Murray  and  Epstein,  1976)  .  Damage  to  such  structures 
may  result  in  a  deterioration  of  the  containment  function  and/or 
danger  to  the  systems  inside  the  structure.  The  limit  states  which 
may  lead  to  such  problems  have  been  identified  by  Murray  and  Epstein 
(1976)  as 

1.  The  pressure  which  causes  cracking  over  a  significant 
portion  of  the  internal  surface. 

2.  The  pressure  which  causes  yielding  of  the  reinforcing 
steel . 

3.  The  pressure  at  which  cracks  may  penetrate  the  structure. 

4.  The  pressure  at  which  the  prestressing  tendons  may  yield, 
thus  seriously  damaging  the  ability  of  the  structure  to 
reseal  itself  upon  relief  of  pressure. 

5.  The  pressure  at  which  spalling  of  concrete  may  lead  to 
debris  inside  the  structure. 

6.  The  pressure  causing  rupture  of  the  reinforcing  steel  or 
prestressing  tendons. 

7.  The  pressure  which  may  initiate  a  structural  mechanism. 

8.  The  pressure  at  which  relative  displacements  may  damage 
the  systems  inside  the  structure. 
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These  limit  states  may  or  may  not  occur  in  the  order  they 
are  listed.  In  Table  6.3  some  of  these  limit  states  are  identified. 
Ultimate  failure  of  this  type  of  structure  may  be  characterized  by 
rupture  of  the  prestressing  tendons.  This  stage  has  not  been  reached 
in  FEPARCS5  analysis.  It  is  possible  however  to  predict  rupture  of 
the  prestressing  tendons  by  studying  the  strain  profiles  and  the  rate 
at  which  strain  increases  with  respect  to  pressure  in  the  late  stages 
of  the  analysis.  The  strain  profiles  at  140.0  psi  for  the  meridional 
and  circumferential  prestressing  tendons  are  shown  in  Figs.  6.22  and 
6.23  respectively.  Fig.  6.22  indicates  that  rupture  of  the  meridional 
prestressing  tendons  may  occur  at  the  top  of  the  dome.  Assuming  that 
the  maximum  strain  of  the  dome  tendons  is  .05,  the  based  on  the 
strain  increment  between  137.5  psi  and  140.0  psi  of  internal  pressure, 
rupture  may  be  expected  at  181.4  psi.  Fig.  6.23  shows  two  possible 
points  at  which  rupture  of  the  circumferential  tendons  may  occur; 
the  top  of  the  dome  and  the  middle  of  the  cylinder.  The  latter, 
however,  possesses  higher  strain  rate  loading  to  failure  at  183.5  psi. 

A  discussion  with  Dr.  S.  Simmonds  has  revealed  that  actual  failure  of 
the  test  structure  took  place  in  the  form  of  slippage  of  two  nonadjacent 
circumferential  tendons  leading  to  a  partial  transfer  of  load  to  the 
middle  tendon  which  then  ruptured  ending  the  experiment  at  159.9  psi. 
This  is  consistent  with  the  findings  of  FEPARCS5 ,  but  it  shows  the 
need  to  establish  slippage  criteria  in  programs  used  for  analysis  of 
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Ultimate  Strength 
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(a)  Material  No.  1  (Normal  Concrete) 
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(b)  Material  No.  2  (Gunite) 


Table  6.1  Concrete  Material  Properties  of  Test  Structure 
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Load 

Step 

Description 

CPU 

(Sec) 

No.  of 

Iterates 

X 

r 

URF** 

NS* 

0 

Problem  Preparation 

2.05 

— 

— 

— 

— 

— 

1 

Prestressing 

37.60 

6 

0.001 

0.005 

0.5 

10 

2 

(G  +  F  +  H)*** 

60.88 

10 

0.001 

0.003 

0.5 

10 

3 

20.0  psi 

48.00 

8 

0.001 

0.003 

0.5 

10 

4 

40.0  psi 

46.65 

8 

0.001 

0.001 

0.5 

10 

5 

50.0  psi 

73.08 

17 

0.001 

0.001 

0.5 

10 

6 

60.0  psi 

65.34 

17 

0.001 

0.003 

1.0 

10 

7 

70.0  psi 

88.05 

23 

0.001 

0.005 

1.0 

10 

8 

75.0  psi 

50.79 

21 

0.001 

0.005 

1.0 

5 

9 

80.0  psi 

48.34 

20 

0.001 

0.007 

1.0 

5 

10 

85.0  psi 

49.37 

21 

0.001 

0.007 

1.0 

5 

11 

90.0  psi 

47.65 

20 

0.001 

0.008 

1.0 

5 

12 

95.0  psi 

43.52 

18 

0.001 

0.008 

1.0 

5 

13 

100.0  psi 

46.12 

19 

0.001 

0.008 

1.0 

5 

14 

105.0  psi 

42.98 

18 

0.001 

0.008 

1.0 

5 

15 

110.0  psi 

47.47 

20 

0.001 

0.008 

1.0 

5 

16 

115.0  psi 

59.61 

26 

0.001 

0.008 

1.0 

5 

17 

117.3  psi 

55.04 

24 

0.001 

0.008 

1.0 

5 

18 

120.0  psi 

56.98 

25 

0.001 

0.009 

1.0 

5 

19 

122.5  psi 

79.14 

33 

0.001 

0.010 

1.0 

5 

20 

125.0  psi 

71.16 

30 

0.001 

0.010 

1.0 

5 

21 

127.5  psi 

79.24 

33 

0.001 

0.011 

1.0 

5 

22 

130.0  psi 

91.20 

39 

0.001 

0.012 

1.0 

5 

23 

132.5  psi 

107.90 

48 

0.001 

0.013 

1.0 

5 

24 

135.0  psi 

110.09 

49 

0.001 

0.015 

1.0 

5 

25 

137.5  psi 

86.06 

38 

0.001 

0.017 

1.0 

5 

26 

140.0  psi 

144.29 

64 

0.001 

0.019 

1.0 

5 

Table  6.2  Load  Steps  of  FEPARCS5  Analysis  of  Finite 
Element  Model  of  Test  Structure 

***  G  :  gravity  load 

F  :  friction  losses 
H  :  hydrostatic  pressure 
**  URF:  under- relaxation  factor 
*  NS:  number  of  subincrements 
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Limit  State  Description 

Pressure 

Location 

extensive  vertical  cracks 

60.0 

cylinder  wall 

extensive  horizontal  cracks 

70.0 

cylinder  and  dome 

ring  beam  cracks 

50.0 

horizontal  cracks  on 

inside  surface 

yielding  of  dome  circumferential 

6  mm  rebars  (72.5  ksi) 

90.0 

top  of  dome 

yielding  of  wall  circumferential 

3  mm  rebars  (48.0  ksi) 

117.5 

from  2'  to  8.5'  above 
base 

yielding  of  dome  meridional  6  mm 
rebars  (72.5  ksi) 

90.0 

top  of  dome 

yielding  of  wall  meridional  6  mm 
rebars  (48.0  ksi) 

110.0 

inside  surface  immedi¬ 
ately  below  ring  beam 

yielding  of  wall  circumferential 
prestressing  tendons  (240.0  ksi) 

130.0 

from  3'  to  7.0'  above 

base 

yielding  of  dome  circumferential 
prestressing  tendons  (220.0  ksi) 

110.0 

top  of  dome 

yielding  of  wall  meridional 
prestressing  tendons  (240.0  ksi) 

— 

yielding  of  dome  meridional 
prestressing  tendons  (220.0  ksi) 

120.0 

top  of  dome 

spalling*of  concrete 

115.0 

inside  surface  of  dome 
near  ring  beam 

Table  6.3  Pressures  Corresponding  to  Limit  States  of  the 
Finite  Element  Model 


* 


spalling  is  defined  as  cracking  parallel  to  the  surface  and  may 
cause  concrete  to  separate  from  the  structure 
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T-Dome  Prestressing  5/e"  Dia.  -  7  W're  Strand 


Fig.  6.1  Elevation  of  Test  Structure 
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Fig.  6.2  Prestressing  and  Rebars  for  the  Dome  of  the 
Test  Structure 
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Fig.  6.3  Detail  of  Ring  Beam  of  Test  Structure 


Fig.  6.4  Detail  of  Cylinder  Wall  Base  Connection  of 
Test  Structure 
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Fig.  6.5  Finite  Element  Mesh  of  Test  Structure 
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(a)  Top  of  Dome 


(b)  Cylinder  Base  (c)  Cylinder  Base 


Fig.  6.6  Boundary  Condition  Simulation  for  Finite 
Element  Model  of  Test  Structure 
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Fig.  6.7  Detail  of  Finite  Element  Model  of  Ring  Beam  Area 
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Fig.  6.8  Reinforcing  and  Prestressing  Layers  of  Finite  Element  Model  of 
Test  Structure 
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Finite  Element  Mesh 


STRSS  (KSI)  STRESS  (KSI) 

0.0  25. 0  50.0  75.0  100.0  0.0  26.0  50.0  75. 
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STRAIN 

(a)  #3  Bars 


0.00  0.02  0.04  0.06  0.08 

STRAIN 

(b)  6  mm  Bars 


Fig.  6.10  Stress  Strain  Curves  for  Prestressing  Bars 
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STRAIN 

(a)  0.5"  (j)  Strand 


STRAIN 

(b)  0.62"  (j)  Strand 


Fig.  6.11  Stress  Strain  Curves  for  Prestressing  Tendons 
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DEFLECTION  (IN) 


Fig.  6.12  Comparison  of  Deflection  Profiles  at  120.0  psi 
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Fig.  6.13  Comparison  of  Internal  Pressure  versus  Maximum  Deflection  of 
Wall  and  Dome 
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Fig.  6.14  Comparison  of  Meridional  Strains  in  Wall  at  28.0  inches  above  Base 
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Fig.  6.15  Comparison  of  Circumferential  Strains  in  Wall  at  28.0  inches 
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Fig.  6.16  Comparison  of  Meridional  Strains  in  Wall 
at  65.0  inches  above  Base 
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.17  Comparison  of  Circumferential  Strains  in  Wall  at 
65.0  inches  above  Base 
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Fig.  6.18  Comparison  of  Meridional  Strains  Near  Top 
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Fig.  6.19  Comparison  of  Circumferential  Strains  Near 
Top  of  Dome 
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Fig.  6.20  Progress  of  Vertical  Cracking 
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Fig.  6.21  Progress  of  Horizontal  Cracking 
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STRAIN 


Fig.  6.22  Profile  of  Strain  in  Meridional 

Prestressing  Tendons  at  140.0  psi 
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STRAIN 

Fig.  6.23  Profile  of  Strain  in  Circumferential 
Prestressing  Tendons  at  140.0  psi 


CHAPTER  SEVEN 


SUMMARY  AND  CONCLUSIONS 


A  nonlinear  elastic  (hypoelastic)  constitutive  relation  has 
been  proposed  for  analysis  of  plane  and  axi symmetric  reinforced  and/or 
prestressed  concrete  structures.  This  constitutive  relation  is  based 
on  the  equivalent  uniaxial  strain  concept  which  has  been  introduced 
by  Darwin  and  Pecknold  (1974,  1977a  and  1977b) .  A  new  characteriza¬ 
tion  of  Poisson's  ratio  has  been  introduced  and  post  failure  conditions 
have  been  imposed  on  the  failure  surface  such  that  the  determination 
of  the  stress-equivalent  uniaxial  strain  relation  parameters  may  be 
more  realistic. 

Special  isoparametric  finite  elements  have  been  developed 
for  representation  of  reinforcing  bars  and  prestressing  tendons  in 
the  meridional  and  circumferential  directions.  These  elements 
together  with  the  well  known  isoparametric  axisymmetric  serendipity 
finite  element  family  (Zienkiewicz ,  1971)  have  been  incorporated  in 
a  finite  element  program  (FEPARCS5)  for  nonlinear  analysis  of  plane 
or  axisymmetric  reinforced  and/or  prestressed  concrete  structures. 

The  proposed  constitutive  relation  has  been  incorporated  in  program 
FEPARCS5  as  well  as  a  number  of  features  such  as  a  post-tensioning 
simulation  procedure. 

Program  FEPARCS5  has  been  used  to  analyse  two  prestressed 
wall  segments  under  biaxial  tension  (Chitnuyanondh ,  et  al.,  1979) 
with  results  comparing  favorably  with  the  test  results  and  with  B0S0R5 
results.  In  analysing  these  segments  the  numerical  instability  of  the 
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original  constitutive  relation  (Elwi  and  Murray,  1979)  has  been 
eliminated  by  the  addition  of  controls  on  the  post  failure  behavior. 
Convergence  has  been  good  up  to  the  yielding  of  the  prestressing  tendons. 
In  these  analyses  the  tangential  stiffness  method  has  been  used  as  an 
iteration  scheme  with  re-evaluation  of  the  stiffness  matrix  every 
three  iterations. 

Finally,  program  FEPARCS5  has  been  used  to  analyse  a  finite 
element  model  of  the  test  structure  erected  and  tested  under  internal 
pressure  at  the  I.F.  Morrison  Laboratory  in  1978  by  J.G.  MacGregor, 

S.  Simmonds,  D.W.  Murray  and  S.  Rizkalla.  The  results  of  the  FE PARCS 5 
analysis  have  been  compared  with  the  test  results  and  with  the  results 
of  an  elastic  plastic  analysis  of  B0S0R5  (Murray,  et  al.,  1978).  The 
comparisons  have  been  based  on  progress  reports  since  the  final 
reports  on  the  test  and  the  B0S0R5  analysis  commissioned  by  the 
Atomic  Energy  Control  Board  of  Canada  have  not  been  published  at 
the  time  of  writing. 

The  constitutive  relationship  developed  herein  appears  to 
be  adequate  for  the  nonlinear  analysis  of  prestressed  concrete 
structures  of  the  type  under  consideration.  On  a  material  level  the 
relationship  has  shown  reasonable  agreement  with  the  two-dimensional 
test  results  of  Kupfer,  Hilsdorf,  and  Riisch  (1969)  ,  as  indicated  in 
Figs.  2.9  to  2.11,  and  with  the  three-dimensional  test  results  of 
Shickert  and  Winkler  (1977),  as  indicated  in  Figs.  2.13  to  2.15. 

When  combined  with  the  strain  softening  characteristics  in  tension, 
shown  in  Fig.  2.2,  it  allows  a  reasonable  simulation  of  the  nonlinear 
response  of  cracked  concrete  under  biaxial  stress  conditions,  as 
indicated  by  the  comparisons  in  Figs.  5.7  to  5.9. 
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However,  a  number  of  problems  with  the  constitutive  model 
require  further  investigation.  More  (numerical)  testing  is  required 
to  see  if  the  softening  branch  is  significantly  affected  by  the 
reinforcing  ratio.  The  tendency  of  a  multiaxial  stress  state  to 
degenerate  to  a  uniaxial  state  under  iteration,  which  has  been  curbed 
herein  by  adopting  the  failure  mode  constraints  described  in  Sect. 
2.3.8  require  further  study.  The  model  should  also  be  compared 
against  material  tests  in  which  the  principal  stress  directions  do 
not  remain  constant. 

The  FEPARCS5  program,  which  has  been  developed  in  associ¬ 
ation  with  this  work,  has  been  shown  to  be  capable  of  reasonably 
simulating  the  nonlinear  behavior  of  a  prestressed  concrete  secondary 
containment  structure  and  hence  to  be  suitable  for  the  analysis  of 
thin-walled  structures  of  this  type.  Prestressing  effects  and  non¬ 
linear  (time-independent)  response  are  treated  without  major  problems. 
The  primary  difficulties  are  the  determination  of  convergence 
criteria  and  the  deterioration  of  the  convergence  characteristics  when 
a  substantial  portion  of  the  concrete  in  the  structure  has  reached 
the  strain-softening  range.  The  latter  difficulty  has  been  dealt 
with  by  using  the  tangent  stiffness  matrix  at  20  psi  as  the  matrix 
for  iteration  of  unbalanced  forces  at  higher  pressures.  However, 
the  entire  area  of  convergence  criteria  and  solution  procedures  for 
structures  containing  strain-softening  materials  is  one  which 
appears  to  require  further  study. 
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APPENDIX  A 


PARAMETERS  OF  ULTIMATE  SURFACES 

As  mentioned  in  Section  2.3.7  the  values  of  the  coefficients 
in  Eqs .  2.43  are  chosen  such  that  the  principal  radii  of  the  ultimate 
surface  and  r2  pass  through  a  set  of  control  points.  for 
C2  =  ^3  <?i  ,  then  0=o°  (Fig.  2.6b)  and  the  points  on  the  surface 

trace  the  variation  of  rj  as  aa  varies.  This  curve  must  pass  through 
the  uniaxial  tensile  strength  point  where  G2  =  a3  =  o  and  Gj  =  f^ 

(and  hence  is  called  the  extension  branch) .  Ic  mu^t  also  pass  through 
the  biaxial  compression  point  where  O 2  =  0 3  =  f  and  Gj  =  o.  The 
nondimensional  values  of  the  tensile  and  biaxial  compression  strengths 
may  be  defined  as 

a  =  f  /f  (A. la) 

t  t  cu 

a  =  f  ,/f  (A. lb) 

c  cb  cu 

One  additional  point  on  the  extension  curve  is  required  to  define  the 
parabolic  variation  of  rx  for  Eq.  2.43a.  An  arbitrary  high  compression 
point  with  G  =  E,  "may  be  selected  from  experimental  data  to  serve 

3 

this  purpose.  Let  the  value  of  T  at  this  point  be  px  as  illustrated 
in  Fig.  2.6c. 

To  define  the  variation  of  r2  (0  =  60°  on  Fig.  2.6b)  it  is 
required  that  Gj  =  a3  _>  G2 .  This  curve  must  pass  through  the  uni¬ 
axial  compression  point  Gj  =  G 3  =  o  and  G2  =  -f  (an<3  hence  is  called 
the  compression  branch) .  It  is  also  required  to  intersect  the  hydro¬ 
static  axis  on  the  rendulic  plane  at  the  same  location  as  the 
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extension  curve  (G  =  £  ) ,  and  to  pass  through  the  experimental 
point  (£,  p2)  where  p2  is  determined  in  the  same  manner  as  px. 

The  values  of  the  coefficients  in  Eqs.  2.43  can  be  deter¬ 
mined  from  the  5  control  points  described  above  and  the  additional 

condition  of  intersection  at  £  as  follows 

o 


a 


2 


a 


l 


a 

o 


b2 


/l  .2  £  (a  -  a  )  -  /l .2  a  a  +  p,  (2a  +  a  ) 

_  t  c _ t  c  1  c  t  #  9 

(2a  +  a)(3£-2a  )(3£  +  a  ) 

c  t  c  t 

(A. 2) 


1 

3 


(2a  -  a) 
c  t 


/l  .2  (a^  -  a  ) 

t  c 

a,,  +  - — - : — 

^  (2a  +a  ) 

c  t 


(A. 3) 


2 

—  a 


4  2 

—  a  a, 
9  c  - 


a 

c 


(A. 4) 


-(a,  +  /  a?  -  4a  a„)/2a,  (A. 5) 

i  1  o  ^  *- 

p2  (£  +1/3)  -  /2/15  (£  +  £) 

(£  +  £  )  (3£-l)  (3£  +1)  *  9  (A'6) 

o  o 


bx  =  (£  +  l/3)b2  +  (/l72  -  3p2)/(3£-l) 


(A. 7) 


b  =  -£  b  -  £2  b  (A. 8) 

o  o  1  ■  o 

The  ultimate  surface  for  stresses  is  therefore  completely 
defined  and  serves  to  evaluate  the  G.  's  required  in  Eqs.  2.21  to 

1C 

2.26. 

However  Eqs.  2.21  to  2.34  also  require  the  evaluation  of  the 
equivalent  uniaxial  strains  associated  with  the  ultimate  strength 

points  G.  as  mentioned  in  Section  2.3.7.  This  is  done  through  the 

1C 
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definition  of  an  equivalent  uniaxial  strain  surface  analogous  to  the 
stress  ultimate  surface.  Control  points  for  strain  can  be  defined 
analogous  to  those  described  for  stresses  and  Eqs.  A.l  to  A. 8  can 
be  used  directly  with  the  appropriate  change  of  variables.  It  should 
be  noted,  however,  that  uniaxial  strains  are  available  from  tests 
only  for  control  points  corresponding  to  f^  and  f^.  Equivalent 
uniaxial  strains  at  the  other  three  control  points  are  fictitious 
strains  that  cannot  be  observed  directly  (Elwi  and  Murray,  1979) . 
These  points  have  been  determined,  herein,  by  trial  and  error  until 
reasonable  strain  correspondence  was  obtained. 
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APPENDIX  B 


ISOPARAMETRIC  SHAPE  FUNCTIONS 


The  shape  functions  for  the  linear,  quadratic  and  cubic 
rectangular  isoparametric  elements  used  in  program  FEPARCS5  are 
defined  in  terms  of  the  normalized  coordinates  y  and  E,  as  follows 
(Zienkiewicz ,  1971). 

(i)  Linear  Elements: 


(l  +  ^i)  (l  +  yyi)/4  ,  i=l ,2,3,4 

(B.l) 

(ii) 

Quadratic  Elements: 

li 

•H 

(l  +  ^i)  (l  +  yyj  (^q  +  yyi-D/4  ,  i=l ,3,5,7 

(B .  2) 

♦i  - 

(1-  C2)  (1  +  yy±)/2  ,  i=2  ,6 

(B.  3) 

ii 

•H 

(1  +  55-.)  (1  -  y2)  /2  ,  i=4 ,8 

(B .  4) 

(iii) 

Cubic 

Elements : 

li 

•H 

(1  +  ££  )  (1  +yyj  [9  (E,z  +y2)  -  10]/32  ,  i=l,4,7, 

10  (B. 

d>.  =  9  (l  +  ££.)  (1  -y2)  (1  +  9yy.  )/32  ,  1=5,6,11,12  (B.6) 

i  i  i 

d>.  =  9(l+yy.)  (l-^2)  (l+9^.)/32  ,  i=2,3,8,9  (B.7) 

l  l  i 

where  i  refers  to  the  node  numbers  appearing  in  Fig.  3.2,  and  E,^  and 
y.  are  the  nodal  normalized  coordinates. 
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APPENDIX  C 


MODELLING  OF  DOME  PRESTRESSING  TENDONS 


The  prestressing  tendons  in  the  dome  of  the  test  structure 
are  arranged  in  two  orthogonal  geodetic  meshes  as  shown  in  Fig.  6.2. 

In  order  to  model  those  tendons  in  terms  of  stiffness  as  well  as 
forces  and  stresses  for  use  in  an  axisymmetric  finite  element  program, 
the  geodetic  meshes  must  be  transformed  into  an  equivalent  meridional 
and  circumferential  mesh.  This  task  is  performed  using  vector 
analysis . 

Let  subscript  £  denote  tendons  which  lie  in  planes  passing 
through  the  y  axis  and  let  subscript  ]i  denote  those  tendons  which 
lie  in  planes  passing  through  the  x-axis  as  shown  in  Fig.  C.l. 

Let  0^  denote  the  angle  plane  describes  with  plane  y-z  and  let  0 
denote  the  angle  that  plane  y  describes  with  plane  x-z.  The  coordin¬ 
ates  of  the  point  of  intersection  of  two  tendons  is  written  as 


X  = 

z  tan 

95 

(C. la) 

y  = 

z  tan 

0 

y 

(C. lb) 

z  = 

Ra 

(C. lc) 

R  is  the 

radius 

of  the  dome  and 

a  =  1.0/(tan2  0>.  +  tan2  0  +  1.0)  l//2  (C.2) 

y 

-)■  ,  , 

Let  P  be  the  unit  vector  along  the  position  vector  of  the 

point  of  intersection  of  two  tendons. 
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P  =  CL  tan  0*.  i  +  a  tan  9  i  +  a  k 

£  y 


(C.  3) 


where  i,  j  and  k  are  the  base  vectors  of  system  x-y-z.  The  unit 
normals  to  planes  E,  and  ]i  are  written  respectively  as 


0, 


cos  0£.  i  +  sin  0^  k 


(C . 4a) 


-> 

0 


y 


cos  0  j  +  sin  0  k 

y  y 


(C.4b) 


The  vectors  tangent  to  tendons  E,  and  y  at  the  point  of  intersection 
are  written  as 


-*■ 


Qr  =  0  x  P 


(C .  5a) 


Q  =  0  x  P 

y  y 


(C .  5b) 


where  x  a  cross  product.  Using  the  definitions  of  Eqs.  C.3  and  C.4, 
Eqs.  C.5a  and  C.5b  can  be  written  as 


Qr  =  (3r  sin  tan  0  )  i  -  (3c-  cos  0r  (tan2  0r  +  1.0)  )j 

9  9  9  y  9  9  9 


+  (3?-  cos  0^  tan  0  )  k 

£  £  y 


(C . 6a) 


-> 

Q. 


-> 


y 


-  (3  cos  0  (tan2  0  +  1.0))  i  +  (3  sin  0  tan  0_)  j 


4-  (3  cos  3  tan  0r)  k 

y  y  £ 


(C.6b) 


where  3 r  and  3  are  defined  as 

9  y 


3  =  1.0/(tan2  0^  +  1.0/cos2  0^)  1//2 


(C.7a) 


3  =  1 . 0/ (tan2  0r  +  1.0/cos2  Q)1^2 

y  9  y 


(C. 7b) 
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The  unit  vector  in  the  circumferential  direction  at  the  point  of 
intersection  of  two  tendons  is  written  as 

C  =  k  x  p  (C.8) 

substituting  for  P  from  Eq.  C.3,  Eq.  C.8  can  be  written  as 


c  = 

Y  tan  0  i  -  y  tan  0,_  j  (C.9) 

h1  S 

where , 


y  = 

1 . 0/ (tan2  0r  +  tan2  0  )1/2  (C.10) 

5  y 

-> 

The  unit  vector  M  lying  along  the  meridian  at  the  point  of  inter¬ 
section  is  written  as 


-> 

M  = 

Cxp  (C.ll) 

Substituting 

for  $  and  P  using  Eqs .  C.3  and  C.9,  Eq.  C.ll  is  written  as 

-*■ 

M  = 

6  tan  0r  i  +  6  tan  0  j  -  5 (tan2  0*.  +  tan2  0  )  k 

£  y  J  K  y 

(C. 12) 

where , 


6  = 

1.0/(tan2  0j-  +  tan2  0  +  (tan2  0-  +  tan2  0  )2)1//2 

s>  y  9  y 

(C. 13) 

The  contributions  of  tendons  £  and  y  to  the  circumferential  direction 
can  be  written  respectively  as 


w  r 
ct, 

=  Q  •  C  (C . 14a) 

€ 

w 

cy 

=  Q  •  C  (C . 14b) 

*y 
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The  contribution  of  tendons  E,  and  y  to  meridional  direction  can  be 
written  respectively  as 

Wm£  =  \  *  M  (C. 15a) 

w  =  Q  •  M  (C. 15b) 

my  t, 

Eqs.  C.14a  to  C.15b  can  be  written  as 


w  r 

=  Y^ 

(sin 

ec 

tan 

2  e 

y 

+ 

tan 

CD 

(cos 

CD 

nnr 

(c. 

.16a) 

w 

=  Y3.. 

(sin 

0 

tan 

2  er 

+ 

tan 

6 

(cos 

e  ) 

(C. 

,16b) 

cy 

y 

y 

£ 

y 

y 

=  a 2  6 

6S 

tan 

% 

cos 

05 

(c. 

,16c) 

w 

=  a2  6 

e , 

tan 

6r 

cos 

e 

(C. 

.  16d) 

my 

y 

£ 

y 

Eqs .  C.16a  to  C.16d  describe  weight  functions  to  be  applied  at  the 
point  of  intersection  of  two  tendons.  The  required  form  of  contri¬ 
bution  must  be  an  average  per  unit  width  for  the  meridional  direction 
and  an  average  per  unit  length  for  the  circumferential  direction.  For 
this  purpose  Eq.  C.16d  has  been  plotted  for  the  6  tendons  appearing 
in  Fig.  6.2  for  mesh  y  as  shown  in  Fig.  C.2.  The  average  meridional 
contribution  was  then  obtained  using  the  formula 

w  =  (w.+2(Zw  +...+w  )  )  /  (7T  R/  2 )  ( C .  1 7 ) 

ma  m1  m2  m6 

This  average  was  then  plotted  on  the  same  figure.  Since  both 
stiffness  and  prestressing  force  as  described  in  Sections  3.3.2  and 
4.4.2  are  directly  related  to  the  area  of  a  meridional  tendon,  the 
average  weight  function  obtained  by  Eq.  C.17  can  be  applied  directly 


to  a  tendon  area. 


. 
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To  obtain  the  contribution  to  the  circumferential  direction 
Eq.  C.16b  has  been  plotted  for  the  6  tendons  appearing  in  Fig.  6.2 
for  mesh  y  as  shown  in  Fig.  C.3.  The  average  meridional  contribution 
was  then  obtained  for  each  ten  inch  arc  length  using  the  formula 

w  =  (1.0  +  £ (w  +  •••  w  )  ) /n  (C.18) 

ca  c  2  cn 

where  n  is  the  number  of  tendons  appearing  in  the  band  described 
by  a  particular  10  inch  arc  length.  These  contributions  were  then 
plotted  on  the  same  figure.  As  with  the  meridional  contributions, 
the  circumferential  weight  functions  can  be  applied  directly  to  the 
area  of  a  tendon  to  simulate  both  stiffness  and  prestressing  force 
formulated  in  Sections  3.3.3  and  4.4.2  respectively. 
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Fig.  C.l  Vector  Representation  of  the  Dome 
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Fig.  C . 2  Meridional  Contribution  of  the  Dome  Prestressing  Mesh 


155 


Fig.  C.3  Circumferential  Contribution  of  the  Dome  Prestressing  Mesh 


